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Abstract 


Recent  studies  of  liquid  films  driven  by  competing  forces  due  to  surface  ten¬ 
sion  gradients  and  gravity  reveal  that  undercompressive  traveling  waves  play  an 
important  role  in  the  dynamics  when  the  competing  forces  are  comparable.  In  this 
paper  we  provide  a  theoretical  framework  for  assessing  the  spectral  stability  of 
compressive  and  undercompressive  traveling  waves  in  thin  film  models.  Associated 
with  the  linear  stability  problem  is  an  Evans  function  which  vanishes  precisely  at 
eigenvalues  of  the  linearized  operator.  The  structure  of  an  index  related  to  the 
Evans  function  explains  computational  results  for  stability  of  compressive  waves. 
A  new  formula  for  the  index  in  the  undercompressive  case  yields  results  consistent 
with  stability. 

In  considering  stability  of  undercompressive  waves  to  transverse  perturbations, 
there  is  an  apparent  inconsistency  between  long-wave  asymptotics  of  the  largest 
eigenvalue  and  its  actual  behavior.  We  show  that  this  paradox  is  due  to  the  unusual 
structure  of  the  eigenfunctions  and  we  construct  a  revised  long- wave  asymptotics. 
We  conclude  with  numerical  computations  of  the  largest  eigenvalue,  comparisons 
with  the  asymptotic  results,  and  several  open  problems  associated  with  our  find¬ 
ings. 
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1  Introduction 


Driven  films  exhibit  a  variety  of  complicated  dynamics  ranging  from  rivulets  and  saw¬ 
tooth  patterns  in  gravity  driven  flows  [JdB92,  SV85]  to  patterns  in  spin  coating  [FH94] 
and  surfactant  driven  films  [TWS89].  A  theoretical  framework  for  these  problems  is 
provided  by  a  lubrication  approximation  of  the  Navier-Stokes  equations  [MW99,  Gre78]. 
This  yields  a  single  partial  differential  equation  for  the  film  thickness  as  a  function  of 
position  on  the  solid  substrate  and  time. 

For  directionally  driven  films,  the  driving  force  enters  into  the  lubrication  approxi¬ 
mation  as  the  flux  /  in  a  scalar  hyperbolic  conservation  law  ut  +  (f(u))x  —  0.  Here  u  >  0 
is  the  film  thickness  and  x  is  the  direction  of  the  driving  force.  For  gravity  driven  flow  on 
an  incline,  the  tangential  component  of  gravity  yields  a  flux  proportional  to  u 3  [Hup82]. 
For  surface-tension  gradient  driven  flows  /  is  proportional  to  u2  [CHTC90].  In  each  of 
these  cases,  the  flux  is  convex.  Consequently,  driven  fronts  in  the  film  correspond  to 
compressive  shock  solutions,  whose  simplest  form  is 


u(x, t ) 


■u_,  if  x  <  st , 

u. |_,  if  x  >  st, 


(1.1) 


in  which  the  shock  speed  s  —  (f(u_)  —  /(i£+))/(ir_ 


m+)  satisfies  the  entropy  condition 


/'(“+)  <  s  <  /'(«-); 


(1.2) 


equivalently,  characteristics  enter  the  shock  on  each  side.  Small  variations  in  height  near 
a  compressive  shock  are  propagated  towards  the  shock  from  both  sides. 

In  practice,  the  discontinuous  fronts  (1.1)  are  smoothed  by  diffusive  effects,  primarily 
through  surface  tension.  In  the  lubrication  approximation,  surface  tension  appears  as  a 
fourth  order  nonlinear  regularization  of  the  conservation  law,  but  there  is  also  second 
order  nonlinear  diffusion  induced  by  the  component  of  gravity  normal  to  the  incline, 
leading  to  the  equation 

ut  +  ( f(u))x  =  -qV  •  (u3VAk)  +  /3V  •  (u3Vu).  (1.3) 

In  this  equation,  7  >  0  and  (3  >  0  are  constants.  The  shock  waves  (1.1)  correspond  to 
smooth  traveling  wave  solutions  of  (1.3);  for  small  (3  >  0,  they  typically  have  oscillatory 
overshoots  and  undershoots  on  either  side  of  the  shock.  Additionally,  the  nonlinearity 
in  the  fourth  order  diffusion  causes  a  single  very  pronounced  overshoot  or  ‘bump’  on  the 
leading  edge  of  the  shock  (see  Fig.  la);  this  structure  is  often  referred  to  as  a  capillary 
ridge  in  experiments.  Capillary  ridges  produced  by  surface  tension  are  well  known  to 
be  linearly  unstable  to  long-wave  perturbations  in  the  transverse  direction  of  the  flow, 
producing  the  well-known  fingering  instability  [CC92,  THSJ89,  KT97,  YC99].  The  effect 
of  larger  (3  is  to  suppress  the  bump  (see  Fig.  lb).  The  disappearance  of  the  bump  (for  (3 
sufficiently  large)  is  accompanied  by  a  transverse  stabilization  of  the  wave  [BB97]. 
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Figure  1:  Compressive  traveling  waves  for  the  flux  /(it)  =  it3,  connecting  it_  =  1  to  u+if 
for  three  values  of  it+.  (a):  (3  —  0,  7  =  1.  All  waves  are  linearly  unstable  to  transverse 
perturbations,  (b):  (3  =  5,  7  =  1.  The  waves  with  it+  =  0.1,0.01  are  linearly  stable  to 
transverse  perturbations. 


Experiments  with  doubly  driven  film  flow,  in  which  gravitational  and  surface  ten¬ 
sion  gradient  stresses  are  competing,  have  uncovered  some  new  phenomena.  While  very 
thin  films  produce  the  characteristic  capillary  ridge  followed  by  a  fingering  instability 
[CHTC90],  thicker  films  produce  a  wider  capillary  ridge  that  continues  to  broaden.  Fur¬ 
thermore,  the  front,  or  leading  edge  of  the  film,  remains  planar;  it  does  not  experience 
fingering  [BMFC98,  KT98]. 

The  lubrication  approximation  for  this  problem  yields  an  equation  of  the  form  (1.3) 
with  a  non-convex  flux 


f(u)  =  u2  -  u3.  (1.4) 

Numerical  experiments  [BMS99,  MB99,  M99]  and  analysis  [BS99]  with  (3  small  show 
that  the  experimental  observations  can  be  explained  by  the  presence  of  undercompres¬ 
sive  shocks.  In  particular,  while  weak  jump  initial  data  give  rise  to  compressive  waves, 
moderately  strong  jumps  can  evolve  into  a  double  shock  structure  in  which  the  lead¬ 
ing  shock  is  undercompressive:  the  speed  s  of  the  shock  violates  the  entropy  condition 
(1.2).  The  trailing  wave  is  compressive  and  travels  at  a  slower  speed.  For  stronger  initial 
jumps,  the  solution  evolves  as  a  rarefaction  wave  connected  to  an  undercompressive  wave 
[MB99];  this  structure  explains  discrepancies  between  earlier  thin  film  experiments  [LL] 
and  analysis  based  on  the  classical  theory  of  conservation  laws,  which  considers  only 
compressive  waves.  The  presence  of  undercompressive  shocks  is  due  to  the  combined 
effects  of  the  non-convex  flux  and  the  fourth  order  diffusion.  Moreover,  numerical  com¬ 
putations  [MB99]  confirm  that  the  undercompressive  leading  wave  is  linearly  stable  to 
transverse  perturbations,  while  the  stability  of  the  compressive  shocks  depends  on  the 
absence  of  a  capillary  ridge. 
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The  numerical  simulations  of  [BMS99]  also  reveal  a  complicated  relationship  between 
initial  condition/far  held  boundary  condition  and  asymptotic  behavior  of  solutions.  More 
specifically,  for  a  range  of  far-held  boundary  conditions,  the  asymptotic  solution  could 
either  be  one  of  a  number  of  compressive  travelling  waves  (see  e.g.  Figure  2),  or  it  could 
approach  the  two-wave  structure  described  above,  with  an  undercompressive  leading 
front.  The  specihc  asymptotic  solution  that  emerges  in  this  case  depends  on  the  shape 
of  the  initial  data  u(x,  0).  For  other  far-held  boundary  conditions,  only  the  two-wave 
structure  is  realized  asymptotically,  irrespective  of  the  shape  of  the  initial  data. 

In  this  paper  we  consider  traveling  wave  solutions  u(x,  y,  t )  =  u(x  —  st )  of  the  general 
equation 


ut  +  Hu)x  =  V  ’  (b(u)Vu)  -  V  •  (c(ii)VAii)  , 


(1.5) 


in  which  we  assume 

(H0)  f,b,ceC2, 

(Hi)  c>  r/  >  0,  6  >  0 

while  keeping  in  mind  the  special  case  studied  in  [BMS99],  corresponding  to: 

f(u )  —  u2  —  u3,  b(u )  =  flu3,  c(u )  =  'yu3,  fl  >  0,7  >  0.  (1-6) 

The  existence  of  traveling  waves  is  itself  an  interesting  problem.  For  convex  fluxes, 
the  existence  of  traveling  waves  can  be  proved  by  a  shooting  argument  as  in  [KH75]  or 
using  topological  methods  involving  the  Conley  index  as  in  [Ren96,  BMS99].  The  case  of 
a  non-convex  hux  such  as  (1.4)  is  more  complicated.  In  particular,  the  phase  portrait  for 
the  resulting  traveling  wave  ODE  can  have  more  than  two  equilibria  and  the  possibility 
of  multiple  heteroclinic  connections,  including  undercompressive  ones.  A  recent  proof 
of  existence  of  undercompressive  waves  for  the  case  (1.5,  1.6)  with  sufficiently  small 
D  —  >  0  was  given  in  [BS99].  Moreover,  non-existence  of  undercompressive  waves 

for  large  D  was  also  established.  Numerical  results  exploring  the  effect  of  varying  D  over 
a  large  range  are  explored  in  [M99]. 

In  this  paper,  we  focus  on  issues  concerning  the  stability  of  travelling  wave  solutions 
of  (1.5).  In  Section  2,  we  consider  stability  to  one-dimensional  perturbations,  and  in 
Sections  3  and  4  we  consider  multidimensional  stability. 

A  stability  theory  based  on  the  Evans  function  [E,  AGJ1,  PW]  has  been  developed  for 
dynamical  systems.  These  ideas  led  naturally  to  the  study  of  one-dimensional  stability 
of  travelling  waves,  including  systems  of  conservation  laws  with  second  order  diffusion 
[GZ,  ZH]  and  for  scalar  conservation  laws  with  second  order  diffusion  and  third  order 
dispersion  [D,  HZ. 2,  Z.3].  Undercompressive  waves  arise  in  both  cases.  In  this  paper,  we 
construct  the  Evans  function  for  the  scalar  equation  (1.5),  in  which  undercompressive 
waves  are  generated  from  the  non-convex  flux  and  fourth  order  diffusion. 

In  Section  2,  following  the  general  approach  of  [GZ],  we  End  a  formula  for  an  index 
T  related  to  the  Evans  function.  When  T  is  negative,  the  travelling  wave  is  unstable, 
since  the  presence  of  a  positive  growth  rate  is  predicted.  When  positive,  the  index  is 
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Figure  2:  Compressive  travelling  waves  for  (1.6)  with  (3  —  0  and  7  =  1,  connecting 
ph  0.332051  to  w+  =  0.1.  Waves  numbered  1,3,  and  5  are  numerically  stable  in 
one  space  dimension,  while  waves  2  and  4  are  unstable.  All  waves  are  unstable  to 
multidimensional  perturbations.  For  this  special  value  of  u_  there  are  an  infinite  number 
of  compressive  waves  connecting  to  u+. 
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an  indicator  of  stability.  (The  index  only  predicts  the  parity  of  the  number  of  unstable 
eigenvalues  or  positive  growth  rates;  when  the  index  is  positive,  we  are  assured  only 
that  there  are  an  even  number  of  unstable  eigenvalues,  not  that  there  are  none.)  The 
geometry  underlying  this  index  is  used  in  Section  2  to  explain  the  following  numerically 
observed  feature:  when  there  are  multiple  compressive  waves  for  the  same  upstream 
and  downstream  thicknesses,  the  waves  come  in  pairs,  one  stable  and  one  unstable  to 
one-dimensional  perturbations.  An  interesting  open  problem  is  to  give  a  rigorous  proof 
of  stability.  Such  a  proof  would  require  further  bounds  on  potential  unstable  eigenvalues 
(perhaps  by  a  variant  of  Grillakis’  method,  as  in  [PW,  AGJS,  D,  P]). 

For  undercompressive  waves,  the  index  requires  special  treatment  due  to  the  behavior 
of  eigenfunctions  associated  with  the  linearization  about  the  travelling  wave.  In  this  case 
we  derive  a  new  simple  expression  for  the  index  that  we  calculate  easily  from  pictures 
of  numerical  results.  This  reduction  turns  out  to  be  necessary  for  practical  evaluation 
of  the  Evans  function  index  in  the  case  of  traveling  waves  for  higher  order  problems.  It 
is  an  improvement  over  the  general  theory  developed  in  [D,  GZ]  and  directly  extends  to 
arbitrary  order  systems. 

The  difference  between  the  compressive  and  undercompressive  cases  appears  even 
more  forcefully  in  the  consideration  of  multidimensional  stability.  In  Section  3,  we  outline 
the  framework  for  a  study  of  stability  to  two-dimensional  perturbations.  We  focus  on 
long-wave  perturbations,  and  show  numerically  in  Section  4  that  compressive  waves 
are  unstable,  in  agreement  with  long-wave  asymptotics  developed  in  Section  3  (and 
also  in  [BB97,  KT98]).  However,  when  considering  undercompressive  waves,  there  is  a 
mathematical  puzzle  we  call  the  long-wave  paradox :  a  formal  calculation  [KT98,  BB97]  of 
the  low-frequency  (long-wave)  expansion  of  the  growth  rate  of  fingering  perturbations  as 
a  function  of  spatial  frequency  yields  an  explicit  formula  for  the  leading  order  asymptotics 
of  the  growth  rate  of  perturbations  at  long  wave  lengths.  While  this  asymptotic  formula 
agrees  well  with  computed  eigenvalues  in  the  case  of  compressive  waves,  it  fails  for 
undercompressive  waves.  In  Section  3  we  resolve  the  long-wave  paradox  by  showing 
how,  for  undercompressive  waves,  the  formal  calculation  overlooks  an  interesting  feature 
of  the  linearized  problem.  In  Section  4,  the  resolution  of  the  long-wave  paradox  is  used  to 
interpret  numerical  calculations  of  growth  rates,  demonstrating  that  undercompressive 
waves  are  stable  to  transverse  perturbations. 

In  Section  4,  we  also  identify  a  curious  feature  of  the  growth  rates.  For  parameter 
values  for  which  there  is  an  undercompressive  wave,  there  are  also  multiple  traveling 
waves,  and  moreover,  the  double  shock  structure  collapses;  the  trailing  compressive  wave 
and  the  leading  undercompressive  wave  have  the  same  speed.  The  graphs  of  the  growth 
rates  (as  a  function  of  wave  number)  for  the  (unstable)  compressive  waves  approach  the 
maximum  of  the  growth  rates  for  the  undercompressive  wave  and  the  compressive  trailing 
wave.  We  interpret  this  observation  in  terms  of  the  linearized  equation,  by  observing 
that  the  trajectories  of  the  compressive  traveling  waves  approach  a  combination  of  the 
trajectories  for  the  undercompressive  wave  and  the  trailing  compressive  wave. 


7 


2  The  One  Dimensional  Case 


To  begin  with,  we  consider  solutions  u  —  u(x,t )  of  equation  (1.5)  that  are  independent 
of  y  : 

V't  T  f(^)x  —  (b('U,)'Ux) x  {c(^)'U‘xxx )x  •  (2.1) 

Suppose  there  is  a  traveling  wave  solution,  which  (by  adding  a  linear  term  to  /,  f(y)  — 
f(y)  —  sy ,)  we  may  take  to  be  stationary: 

u(x ,  t)  —  u(x)  (2-2) 

Then  u  —  u(x)  satishes  the  traveling  wave  ODE 

c(u)uin  =  b(u)u'  -  ( f(u )  -  /(«+))  (2.3) 

and  boundary  conditions 


lim  u(x)  —  u±.  (2-4) 

x— *-d=oo 

In  addition  to  the  standing  assumptions  (Ho),  (Hi),  we  assume  standard  nondegeneracy 
conditions  from  [ZH]4  : 

(H2)  a±  :=f'(u±)^  0. 

(H3)  Solutions  of  (2.3)  -  (2.4)  are  (locally  to  u)  unique  up  to  translation. 

Remark  (H2)  holds  if  and  only  if  (rx-t  ,0,0)  are  hyperbolic  rest  points  of  the  first  order 
system  (for  u,u',u")  associated  with  (2.3).  Equivalently,  it  states  that  that  the  wave 
speed  is  non-characteristic  for  the  underlying  conservation  law  ut  +  (f(u))x  —  0. 

Linearizing  (2.1)  about  u(-)  gives 

vt  =  Lv  :=  ~(cvxxx)x  +  (bvx)x  -  {av)x, 
in  which  c  —  c(u(x)),  b  —  b(u(x)),  a  —  f  '(u(x))  —  b'(u)ux  +  c'(u)uxxx. 

Letting  v(x,t)  —  eXtiu(x),  we  obtain  the  eigenvalue  problem  Liu  —  Xw  : 

—  ( cw '") '  +  {bw ') '  —  (aw) '  —  Xw.  (2-5) 

Here,  and  below,  '  denotes  differentiation  with  respect  to  x. 

To  demonstrate  linear  stability  of  the  travelling  wave  u,  we  want  to  show  that  if 
Re X  >  0,  A  7^  0,  then  equation  (2.5),  with  suitable  boundary  conditions  at  ±00,  has  only 
the  trivial  solution  w  =  0.  Our  approach  to  this  question  is  through  the  Evans  function 
D( A)  (defined  below)  of  [AGJS,  GZ].  This  is  an  analytic  function  whose  zeroes  in  the 
right  half  plane  (minus  the  exceptional  point  A  =  0,  which  lies  on  the  boundary  of  the 

4(H3)  is  the  specialization  to  the  scalar  case  of  the  more  general  (H3)  in  [ZH]. 
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essential  spectrum)  correspond  precisely  to  eigenvalues  of  L  .  Moreover,  as  shown  in 
[ZH,  HZ. 2],  stability  under  quite  general  circumstances  is  equivalent  to  the  condition: 

(V)  D(-)  has  precisely  one  zero  on  {Re( A)  >  0},  consisting  of  a  simple  root  at  A  =  0. 

The  Evans  function  itself  is  rather  difficult  to  evaluate;  a  more  readily  computable 
quantity  is  the  stability  index  [J,  PW,  GZ,  Z.2,  Z.3] 

T  :=  sgn  D'( 0)  lim  sgn  D{ A),  (2-6) 

A — >-oo 

where  the  limit  is  along  the  real  axis.  This  quantity,  taking  values  T  =  ±1,  gives  the 
parity  of  the  number  of  zeroes  of  D(-)  in  the  unstable  half-plane  Re A  >  0.  This  follows 
from  the  observations  that  (i)  T  records  the  number  of  crossings  of  D( A)  through  zero  as  A 
moves  along  the  positive  real  axis,  and  (ii)  non-real  zeroes  of  D  appear  in  conjugate  pairs, 
since  L  is  real.  In  particular,  T  =  —  1  implies  instability  of  the  traveling  wave,  whereas 
T  =  1  is  consistent  with  stability.  It  is  this  computable  index,  and  its  interpretation 
in  terms  of  the  phase  portrait  of  the  traveling  wave  ODE,  that  provides  much  of  the 
theoretical  explanation  of  the  numerically  observed  stability  phenomena. 

To  construct  the  Evans  function,  we  explore  solutions  of  (2.5)  as  follows.  Since  the 
equation  is  fourth  order,  for  each  A,  there  will  be  a  four  dimensional  set  of  solutions. 
For  A  e  {Re( A)  >  0}/{0},  we  can  choose  a  basis  for  the  subspace  <S+  of  solutions  that 
approach  zero  as  x  — >  oo  and  a  basis  for  the  subspace  U~  of  solutions  that  approach 
zero  as  x  — >  —  oo.  Then  A  is  an  eigenvalue  if  <S+  and  U~  have  non-trivial  intersection. 
That  is,  we  need  a  condition  for  the  intersection  of  two  linear  two-dimensional  subspaces 
of  functions.  This  condition  is  that  the  Evans  function  D( A),  defined  to  be  the  Wronskian 
of  the  two  pairs  of  basis  functions  of  <S+  and  U~,  should  vanish.  In  the  present  situation, 
the  sign  of  the  Wronskian  is  shown  to  be  independent  of  x.  The  construction  of  the 
Evans  function  extends  analytically  to  A  =  0. 

In  the  construction  of  -D(A),  we  do  not  need  to  include  information  about  whether  the 
underlying  traveling  wave  is  compressive  or  undercompressive.  After  the  Evans  function 
is  defined  and  related  to  stability  of  traveling  waves,  we  then  focus  on  the  compressive  and 
undercompressive  cases  in  turn.  In  sections  2.5  and  2.6,  we  show  that  these  cases  differ 
most  significantly  at  A  =  0,  due  to  the  fact  that  small  disturbances  propagate  through 
an  undercompressive  wave,  rather  than  being  absorbed,  as  they  are  for  a  compressive 
wave.  This  distinction  reveals  itself  in  the  structure  of  the  subspaces  <S+  and  U~ ,  but  it 
can  already  be  seen  in  the  behavior  at  x  —  ±oo,  which  we  now  discuss. 


2.1  Asymptotic  eigenvalue  equations. 

As  x  -7-  Too,  behavior  of  (2.5)  is  governed  by  the  asymptotic  constant  coefficient  equa¬ 
tions 


c±w  ""  —  b±w  "  —  a±w '  —  A w. 


(2.7) 
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(a)  A  — »  +oo  (b)  A  =  0,  b±  =  0,  a±  <  0 

Figure  3:  Roots  of  equation  (2.8) 

Here  and  below,  we  use  subscripts  ±  to  indicate  that  quantities  depending  on  u  are 
evaluated  at  u  —  u(dzoo),  respectively.  E.g.,  c+  =  c(u(oo)).  In  particular,  a±  —  f'(u±), 
so  that  the  sign  of  a±  governs  whether  characteristics  for  the  underlying  conservation 
law  point  towards  the  stationary  front,  or  away  from  the  front.  This  distinction  shows 
the  effect  of  the  conservation  law  on  the  direction  of  propagation  of  small  disturbances. 
The  normal  modes  w  :=  x  of  (2.7)  are  determined  by  the  characteristic  equation 

c-t/i4  —  b±p2  +  a±p  +  A  =  0.  (2-8) 

First  note  that  for  A  =  0,  one  of  the  roots  is  /i  =  0.  This  corresponds  to  the  constant 
solution  of  equation  (2.7)  when  A  =  0.  For  A  ^  0,  we  find  that  all  solutions  of  (2.8)  have 
nonzero  real  part. 

Recall  that  we  need  to  distinguish  solutions  w(x)  that  decay  at  x  —  Too,  for  which 
Re(p )  <  0  (termed  stable  in  Lemma  2.1  below)  and  solutions  w (x)  that  decay  at  a;  =  —  oo, 
for  which  Re(p)  >  0  (termed  unstable  in  Lemma  2.1). 

Lemma  2.1  For  A  £  {i?e(A)  >  0}/{0};  (2.7)  has  two  stable  (i.e.  negative  real  part) 
and  two  unstable  roots ,  pi,  p2  <  0  <  /13,  P4  (ordering  by  real  parts).  At  A  =  0,  there 
are  two  cases: 

(a±  >  0)  :  pi  <  p2  -  0  <  p3,  P4- 

(a±  <  0)  :  pi,  p2  <  0  =  pz  <  /i4. 

Proof.  (A  7^  0):  First,  observe  that  (2.8)  has  no  imaginary  roots  p  for  A  £  {Re( A)  > 
0}  \  {0}.  For,  setting  p  —  ki  in  (2.8),  we  have  the  dispersion  relation 

A  =  —a±ik  —  b±k2  —  c±ki,  (2-9) 

giving  Re( A)  <  0  unless  k  —  0.  Thus,  the  number  of  stable /unstable  roots  is  constant. 
Taking  A  — >  +00,  on  the  real  axis,  we  have  p  (— A)1//4  =  | A|1//4( —  l)1/4,  giving  two 
stable  and  two  unstable  roots  (Fig.  3(a)). 

(A  =  0):  At  A  =  0,  we  can  factor  (2.8)  as 

c±p3  —  b±p  +  a±  —  0,  p  —  0.  (2.10) 
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Clearly,  the  first  equation  has  no  imaginary  roots  for  a±  /  0,  so,  we  can  again  count 

/  \l/3 

roots  by  homotopy,  taking  b±  ->•  0  to  find  ji  =  (  —  ~  )  This  gives  two  stable  roots  and 

one  unstable  root  for  a±  <  0,  together  with  /i  =  0,  as  claimed  (Fig.  3(b)).  The  case 
a±  >  0  is  symmetric.  ■ 


Now  write  (2.7)  as  a  first  order  system 


W'  =  A±(A  )W, 


A±  : — 


/ 

V 


0 

0 

0 


—  A/c 


1  0 

0  1 

0  0 

—a/c  b/c 


0  \ 
0 
1 


°/± 


(2.11) 


in  the  phase  variable  W  (w,  iu  ',  iu  ",  w  '")t .  Throughout  this  paper  we  will  identify 
solutions  of  such  fourth  order  scalar  ODEs  with  solutions  of  their  corresponding  first 
order  system. 

For  A  /  0,  each  root  /i*,  is  associated  with  an  eigenfunction  w{x)  —  enx,  a  solution  of 
(2.7).  Then  W  ,iu" ,iu'"Y  is  a  solution  of  (2.11).  Since  the  equation  is  linear, 

these  exponential  solutions  can  be  combined  linearly  to  form  subspaces  of  solutions.  In 
keeping  with  the  strategy  of  distinguishing  exponentially  decaying  from  exponentially 
growing  solutions,  we  separate  the  corresponding  subspaces  of  solutions  S+ ,IA  ,  where 
<S+  is  the  two-dimensional  invariant  subspace  of  solutions  that  decay  as  x  — >  oo,  and 
IA  is  the  two-dimensional  invariant  subspace  of  solutions  that  decay  as  x  — >  —  oo.  The 
content  of  the  following  Corollary  is  that  not  only  are  these  subspaces  two  dimensional 
for  A  ^  0,  Re  (  A)  >  0,  but  they  also  extend  analytically  to  A  =  0  as  two  dimensional 
subspaces,  even  though  for  A  =  0  each  subspace  may  have  a  non-decaying  eigenfunction 
(depending  on  the  sign  of  a±)  as  described  in  Lemma  2.1. 

In  the  case  of  a  compressive  wave,  a_  >  0  and  a+  <  0.  Thus  for  A  =  0  both  S+ 
and  U  are  composed  of  decaying  eigenfunctions  at  their  respective  infinities  in  x.  On 
the  other  hand,  for  an  undercompressive  wave,  a_  and  a+  have  the  same  sign.  If  they 
are  both  negative  then  for  A  =  0,  IA  has  one  nondecaying  eigenfunction  at  x  — >  —  oo 
while  has  two  decaying  modes  at  oo.  This  behavior  carries  over  to  the  non-constant 
coefficient  case  describing  the  full  eigenvalue  problem  (2.5). 

Corollary  2.2  The  s table /unstable  subspaces  ZT1" jU.  associated  with  A-t(-)  are  each 
two-dimensional  on  {Re( A)  >  0}/{0},  and  extend  analytically  in  A  to  {Re A  >  0}. 
More  precisely,  there  exist  bases  {Ob1",  V/ },  {O3-,  V4~ }  for  the  subspaces  S+ ,IA  : 


S*  -  -pm. !  t •  I  V  }, 


(2.12) 


U  —  span{f  b  ,  V4  }. 


(2.13) 


The  bases  can  be  chosen  to  depend  analytically  on  A  and  have  the  symmetry  V/(A)  — 

Wi¬ 


ll 


Proof.  The  dimension  follows  from  Lemma  2.1.  Likewise,  since  groups  /ii,  /i2  and  /i3, 
/i4  remain  spectrally  separated  on  {Re A  >  0},  the  generalized  eigenprojections  onto 
their  associated  subspaces  are  each  analytic,  by  standard  matrix  perturbation  theory 
[K].  The  eigenprojections  are  given  by  the  resolvent  formula  P  :=  fr(A±  —  /ul)~1d/i ,  T  a 
contour  enclosing  /i^,  /if  (resp.  /if,  /i^).  If  T  is  chosen  to  be  symmetric  with  respect  to 
complex  conjugation,  then  P  is  as  well,  by  the  corresponding  property  of  A± .  Choosing 
an  analytic  basis  of  eigenvectors  by  the  construction  of  [K,  pp.  99-102],  we  retain  the 
above  symmetry.  ■ 

2.2  Stable  and  unstable  manifolds 

In  the  previous  subsection,  we  froze  the  coefficients  in  the  ODE  (2.5),  by  setting  u  —  u+ 
and  u  —  u_.  This  allows  us  to  capture  the  behavior  of  solutions  of  the  nonconstant 
coefficient  equation  (2.5)  as  x  — >  ±oo. 

Given  a  solution  <f>(x)  of  (2.5),  we  associate  with  it  the  vector  of  derivatives  $(a;)  = 
(4>(x),  4>"{x),  Given  such  a  vector  $,  the  angle  it  makes  with  a  subspace 

T  of  the  four  dimensional  phase  space  is  the  minimum  of  the  angle  between  $  and  all 
vectors  in  T. 

On  {Re A  >  0},  define  to  be  the  subspace  of  solutions  </>  of  (2.5)  whose  corre¬ 
sponding  vectors  $  approach,  in  angle,  the  subspace  as  x  — >  Too.  Analogously, 
U~  is  defined  to  be  the  subspace  of  solutions  </>  of  (2.5)  whose  corresponding  vectors  $ 
approach,  in  angle,  the  subspace  U~  as  x  — >  —  oo.  The  existence  of  these  subspaces  is 
established  in  [GZ]  or  [K.l,  K.2,  K.3]. 

The  existence  of  these  subspaces,  along  with  analytic  dependence  on  A,  follows  by  the 
gap  lemma  of  [GZ,KS],  or  alternatively  by  earlier  results  of  [J,K.  1— 3] ;  the  full  generality  of 
the  gap  lemma  is  not  needed  here,  since  we  have  spectral  separation  of  the  subspaces  <S+, 
U.  and  the  complementary  A± -invariant  subspaces,  a  helpful  feature  of  the  scalar  case. 
The  necessary  hypotheses  follow  by  analytic  dependence  of  S'"1",  U  ,  the  aforementioned 
spectral  gap,  and  exponential  convergence  of  the  coefficients  a,  6,  and  c  as  x  — >  ±oo. 
Regarding  the  latter  property,  recall  from  Remark  1  that  u±  are  hyperbolic  rest  points 
of  (2.5)  by  virtue  of  (H2),  hence 

\(d/dx)\u-u±)\  <  Ce-aN,  £  —  0,1,2  (2.14) 

for  some  a  >  0,  and  thus 

\a- a±\,\b-b±\,\c- c±\  <Ce~alxl,  x^O.  (2.15) 

Using  a  slight  extension  of  these  results  (see  [ZH]  Lemma3.1,  p.  779,  or  [GZ]  Corollary 
2.4,  p.  807),  we  can  conclude  in  addition  that  there  exist  basis  functions  p±(A,x)  such 
that 


S+  =  span{^,^}, 


(2.16) 
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U  =  span{<p3  ,  (p4  }, 


(2.17) 


which  additionally  have  symmetry  pf(X)  =  ^pf(X)  so  that  are  real- valued  for  A  real. 
This  follows  from  the  properties  noted  above,  plus  the  existence  of  corresponding  basis 
vectors  asserted  in  Corollary  2.2.  Moreover,  the  basis  functions  ft(X,  x)  are  analytic 
in  A  at  A  =  0  and  continuous  elsewhere;  the  exterior  products  (pf  A  pf  and  fz  A  pf  are 
analytic  on  all  of  {Re A  >  0};  and  p(X)  =  p(X)  (see  [GZ]). 

Because  their  associated  eigenvalues  may  coalesce,  <f>f,  </> \  may  not  be  individually 
analytically  continuable  on  {Re A  >  0};  however,  they  are  jointly  continuable  in  the  sense 
that  their  span,  or  alternatively  their  exterior  product,  can  be  analytically  defined. 

2.3  The  Evans  function. 

We  now  define  the  Evans  function,  following  [AGJS,  GZ]  as  the  Wronskian  of  the  basis 
functions  pf,  pif,  p 3  ,  pf  that  characterize  solutions  decaying  at  +00,  at  —00,  or  both. 


/  fl  ft  fz  fi  \ 


D( A)  :=  det 


<pV 


<pi" 
\  ft'" 


///  / 
f4  / 


Properties: 

(PI)  D(-)  is  analytic  on  {Re A  >  0}. 

(P2)  D(-)  is  real  for  real  A. 

(P3)  On  {.ReA  >  0}/{0},  D( A)  =  0  if  and  only  if  A  is  an  eigenvalue  of  L. 


(2.18) 


Only  (P3)  requires  discussion.  This  follows  from  the  characterization  of  eigenfunc¬ 
tions  of  L  as  nontrivial  solutions  of  (L  —  A )w  —  0  lying  in  <S+  0  U~ ,  i.e.  decaying  at  both 
±00.  Evidently,  vanishing  of  the  Wronskian  D(-)  is  equivalent  to  nontrivial  intersection 
of  and  U~,  by  (2.16)  -  (2.17). 

The  meaning  of  the  Evans  function  at  A  =  0  is  less  immediate,  but  equally  important. 
Note  that  A  =  0  is  an  eigenvalue,  with  at  least  the  eigenfunction  ux,  corresponding  to 
translations  of  u.  Moreover,  since  ux  decays  at  ±00,  it  lies  in  both  U~  and  <S+,  implying 
D( 0)  =  0.  Moreover,  this  zero  of  the  Evans  function  does  not  interfere  with  either  linear 
or  nonlinear  stability,  as  shown  in  [HZ. 2,  ZH],  a  result  we  alluded  to  earlier  in  formulating 
condition  (T>),  and  which  we  restate  here  for  clarity. 


Proposition  2.3  Linearized  stability  of  u(-)  as  a  solution  of  (2.1)  is  equivalent  to  the 
Evans  function  condition  (V) .  Moreover,  linearized  stability  implies  nonlinear  stability. 
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(Remark:  [HZ. 2]  concerns  also  more  general,  dispersive-diffusive  equations  such  as 
the  convex  KdV-Burgers  equation,  or  the  nonconvex  modified  KdV-Burgers  equation 
studied  by  [D,  JMS,  W].) 

Having  defined  the  Evans  function,  we  now  begin  to  investigate  whether  or  not  it  has 
a  positive  zero.  To  this  end,  we  use  the  stability  index  T  (see  (2.6))  which  only  requires 
information  about  D( A)  for  sufficiently  large  A  and  the  leading  order  behavior  for  real 
A  near  A  =  0.  The  small  A  behavior  depends  crucially  on  whether  the  traveling  wave  is 
compressive  or  undercompressive.  The  stability  index  will  be  shown  to  be  a  coordinate- 
independent  orientation  of  the  intersection  of  the  stable /unstable  manifolds  of  u+/u_  in 
(2.3).  The  connection  to  stability  is  given  through 

Proposition  2.4  The  parity  of  the  number  of  zeroes  of  £)(•)  in  {i?e(A)  >  0}  is  odd 
(even)  according  to  T  >  (<)0.  In  particular,  T  <  0  implies  instability,  whereas  T  >  0  is 
necessary  for  stability. 

Proof.  Using  D( A)  =  D( A),  we  have  that  complex  roots  of  D(-)  appear  in  conjugate 
pairs.  On  the  other  hand,  the  number  of  real  roots  with  A  >  0  clearly  has  the  parity 
claimed.  The  connection  to  stability  follows  from  Proposition  2.3.  ■ 

2.4  The  Evans  function  as  A  — >•  oo. 

Following  [GZ],  we  now  evaluate  the  Evans  function  in  the  large  |A|  regime.  Let  7r  :  C4  — > 
C4  denote  orthogonal  projection  onto  the  span  of  the  first  two  standard  basis  elements, 
e\  —  (1,0,  0,0)*  and  e2  =  (0, 1,0,0)*.  We  have  the  following  analog  of  Lemma  3.5  in 

[GZ]: 

Lemma  2.5  The  projection  n  is  full  rank  on  the  A±  -invariant  subspaces  S+ ,  U  .  Equiv¬ 
alently,  ifV  —  {v\,  V2,  t>3,  tq)4  and  V  —  (ui,  v2,  £3,  rq)4  £  <S+  (resp.  U  )  are  independent, 

then  det  (  ll  J  7^  0. 

\  v2  v2  J 


Proof.  It  is  sufficient  to  treat  the  case  that  V,  V  are  eigenvectors  of  A±.  If  they  are 
independent  genuine  eigenvectors,  then,  by  the  companion  matrix  structure  of  (2.11),  we 

have  V  =  (l,  fi,  /12,  /i3),  V  =  (1,  fi,  jl2,  (l3),  and  det  ^  ^  ^  ^  =  T  ~  T  +  0. 

The  other  possibility  is  that  y  —  fi  and  V  —  (1,  /i,  /i2,  /i3)4,  V  —  (0, 1,  2/i,  3/i2)4  form 


an  eigenvector  /  generalized  eigenvector  pair.  In  this  case,  det 


V1  V! 
V2  V2 


—  1^0,  also. 


This  lemma  allows  us  to  express  the  sign  of  D(X)  for  large  A  in  terms  of  the  projections 
of  the  basis  vectors  of  ZT1",  U  at  A  =  0.  To  this  end,  let  the  basis  vectors  1 Ip  (j=l,...,4) 
have  components  Vf(.  (k=l,...,4). 
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Proposition  2.6  For  A  >  0  real,  sufficiently  large, 

sgn  D( A)  =  sgn  del  (  ])|  ^  )  det  (  J”  ^  )  U_o-  (2.19) 

Proof.  The  proof  follows  arguments  of  Proposition  2.8  of  [GZ],  or  the  equivalent  Propo¬ 
sition  7.3  of  [ZH]  for  lower  order  systems  of  equations.  The  idea  is  to  use  an  asymptotic 
argument  for  large  A  to  show  that  the  sign  of  the  determinant  in  the  Evans  function  can 
be  expressed  in  terms  of  the  sign  of  the  determinants  in  (2.19)  provided  A  is  sufficiently 
large.  Then  Lemma  2.5,  which  implies  that  the  sign  of  the  determinants  in  (2.19)  are 
independent  of  A  (since  they  are  continuous  and  can  never  vanish)  gives  the  desired 
result. 

First  we  choose  A  sufficiently  large  so  that  the  asymptotic  arguments  below  are  valid. 
We  can  rescale  (2.5)  by  x  —  A1//4m,  w(x)  —  w(x),  a(x)  —  a(x),  b(x)  —  b(x),  c(x )  =  c(x), 
to  obtain 

w""  =  -c_1ru  +  e>(|A“1/4|(H  +  \w'\  +  \w"\  +  \w  ,,;|))  (2.20) 

where  |c'|  =  0( A_1//4).  Applying  Proposition  2.8  of  [GZ],  or  the  equivalent  Proposition 
7.3  of  [ZH],  we  find  that  in  phase  variables,  the  subspaces  S+(x)/U~(x),  defined  by  the 
evaluation  of  the  stable /unstable  subspace  of  (2.20)  at  an  arbitrary  xq,  lies  within  angle 
0( A-1//4)  of  the  stable/unstable  subspaces  of  the  limiting  equations 

w =  —c~1iu,  (2.21) 

with  coefficient  c  held  frozen  at  value  c(x o).  The  subspaces  of  the  limiting  equations 
are  spanned  by  the  stable/unstable  normal  modes  of  (2.21),  readily  calculated  to  be 
(1,  f,  f2,  f3y,  where 

f  =  fj(x  o)  =  (-c(£0)_1)1/4  (2.22) 


are  fourth  roots  of  (— c)_1. 

Converting  back  to  the  original  scaling,  the  subspaces  S+(xq)/U~(xq)  defined  by  the 
evaluation  of  the  stable /unstable  subspace  of  (2.5)  at  x(),  lie  within  angle  0( A-1/4)  of  the 
subspaces  spanned  by  {(1,  ^i,  Mi,  Mi)*,  (1?  ^2,  •  •  •  ,  l4  Y)  and  {(1»  ^3,  t4,  l4)i  (^^4,  &  ^4)} 
respectively,  where 


Vj(x  0)  =  (-A/c(a:o))1/4. 


(2.23) 


Recalling  (see  Section  2.2)  that  S+(xq)/U  (so)  are  real-valued  for  A  real,  and  that 
for  large  A  /i  1,  /i2  and  /i3,  fi^  form  complex  conjugate  pairs,  we  find  therefore  that 
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where  M±(x)  are  specific  nonsingular  real- valued  coefficient  matrices.  It  follows,  there¬ 
fore,  that 


£(A)  =  -(l  +  C>(|A-1/4|))det 


/  1  1 

i  i  > 

1*  1  1*2 

CO 

4^ 

det  | 

'  M+ 

0 

f*i 

0 

V. 

\  i4  ■■■ 

f*4  j 

\x=x0 

(2.25) 


Evaluating  the  Vandermonde  determinant,  and  noting  that  all  quantities  involved 
are  real  and  nonvanishing,  hence  have  sign  independent  of  xq,  we  find  that 


sgn  D( A) 


,  .  •  (I  1 

sgn  det i 

V  ¥l  Ah 

,  .  •  (I  1 

sgn  det i 

\  l*i  1*2 

—  sgn  det 


det  M+  \X=X0  det 


det  M+\x=+oc  det 


v*i  K 

v&  V+ 


det 


1*3  1*4 

1  1 

f*3  1*4 

V£  V£ 
V£  V4~2 


so  long  as  A  is  real  and  sufficiently  large. 

But,  the  right  hand  side  of  (2.26)  is  independent  of  A  by  Lemma  2.5.1 
Remark:  The  matrices  M±  in  the  proof  are  used  only  as  an  intermediate  book¬ 
keeping  device,  and  play  no  role  in  the  final  computation. 


2.5  Evaluation  at  A  =  0  in  the  Lax  case. 

In  this  subsection,  we  restrict  to  the  Lax  case,  for  which 


a_  >  0  >  a+.  (2.27) 

Consulting  Lemma  2.1  and  Section  2.2,  we  find  that  <S+(0),  Z7_(0)  consist  of  all 
solutions  of  (2.5)  that  exponentially  decay  as  x  — >  +oo,  —  oc  respectively.  We  can 
therefore  integrate  (2.5)  with  A  =  0  from  either  +oo  or  — oo  to  obtain 
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cp'"  —  bp'  —  ap.  (2.28) 

Note  that  (2.28)  is  exactly  the  variational  equation  for  the  traveling  wave  ODE  (2.3). 
Indeed,  one  solution  of  (2.28)  is  p  —  ux,  corresponding  to  translation  along  the  profile. 
By  appropriate  change  of  basis,  we  can  arrange  without  loss  of  generality  that 

<p]*~(0,  x)  —  p~l  (0,  x)  =  u(x).  (2.29) 

We  choose  p^  and  p%  to  be  any  independent  solutions  of  (2.28)  decaying  at  ±oo. 

Note  that  the  traveling  wave  ODE  (2.3)  can  be  written  as  a  three  dimensional  au¬ 
tonomous  system 


u 

/ 


Hi, 

U2, 

b(U). 

c(u) 


(2.30) 


Hi 


/(h)  -  /'(h_) 
c(h) 


which  has  equilibria  at  B  —  (h+,0,0)  and  M  —  (h_,0,0).  A  travelling  wave  solution 
of  (2.3)  connecting  h_  to  u+  corresponds  to  a  heteroclinic  orbit  of  (2.30)  connecting 
M  to  B.  For  the  Lax  case,  M  has  a  two  dimensional  unstable  manifold  and  B  has  a 
two  dimensional  stable  manifold.  We  now  see  that  the  subspace  <S+  (U~ )  at  A  =  0  is 
composed  of  functions  p  with  corresponding  vectors  (p,  p',  p"Y  that  span  the  tangent 
space  of  WS(B)  ( WU(M ))  along  h(-).  An  illustration  of  this  structure  for  the  special 
case  of  (1.6)  is  presented  at  the  end  of  this  section. 

The  next  proposition,  following  the  general  approach  introduced  by  Jones  [J]  relates 
the  sign  of  D'( 0)  to  the  orientation  with  which  these  stable /unstable  manifolds  intersect 
in  the  phase  space  M3  of  (2.3). 


Proposition  2.7  D( 0)  =  0,  while 
sgn  D'( 0)  =  sgn 

Proof. 

D{  0)  =  det 


l 

lP3i 

c(O)-1 

det  1 

K 

+' 

V 

K 

+" 

V2 

f3 

^  ^ X 

<P2 

^ X 

III 

\  ux 

(rfy 

ii 

(Vs)'" 

III 

ux 

=  0. 


/ 


Similarly,  a  straightforward  calculation  (as  in  [GZ])  yields 

ux  <^2  ^3  (W4A  -  \ 

D'( 0)  =  det  ' 


+'"  (  -  +  \m  I 

^3  W4A  -<Plx)  ) 


H„ 


(2.31) 


(2.32) 


(2.33) 
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Differentiation  of  (2.5)  with  respect  to  A  shows  that  the  functions  z  —  (pf  ,  <p4  satisfy 
the  same  variational  equation 


(cz'J  =  (bz'Y  -  (az)'  -ux, 


(2.34) 


decaying  at  +oo,  — oo  respectively.  Integrating  from  +oo,  — oo  respectively,  and  sub¬ 
tracting,  we  find  that  z  —  ipf  satishes 

cz"  —  bz  —  az  +  (w_  —  u+),  (2.35) 


an  inhomogeneous  version  of  (2.28).  Using  (2.28),  (2.35)  to  eliminate  the  fourth  row  in 
(2.33),  we  obtain 


D'( 0)  =  c(0)-1  det 


^  ux 

_ n 

Ux 


V  0 


1  / 

V2 

+" 

1 

V2 

Vs 

0 

0 

5  \ 

z’ 

z" 

( U -  -  U+) 


(2.36) 


giving  the  result.  ■ 

Remark:  The  right  hand  side  of  (2.36)  is  of  form  qA,  where  the  first  factor,  7 
measures  orientation  of  the  intersection  of  the  stable  and  unstable  manifolds  of  the 
underlying  traveling  wave  ODE  at  u+,  u_  respectively,  and  the  second  factor,  A  =  [«], 
is  the  Kreiss-Sakamoto-Lopatinski  determinant  arising  in  the  study  of  inviscid  stability. 
This  is  a  special  case  of  a  very  general  relation  pointed  out  in  [ZS]  between  viscous  and 
inviscid  stability. 

Combining  Proposition  2.6,  Proposition  2.7,  we  have 


Corollary  2.8  The  stability  index  T  sgn  D'(0)D(+oo)  satisfies 


T  = 


sgn 


det 


ux 

_ t 

un 


V2 

<P2 


lcc=0 


X 


det 


(p3  ux 
<P3  K 


c  =  —  00  (^—  ^+) 


We  now  present  an  example  from  [BMS99]  and  discuss  how  the  machinery  developed 
here  explains  numerical  observations  of  one  dimensional  stability  of  compressive  waves. 
We  consider  (1.6)  with  (3  =  0  and  7=1.  For  each  u+  there  is  a  range  of  for  which 
multiple  compressive  traveling  waves  appear.  For  a  special  value  of  u_  an  infinite  number 
of  compressive  waves  occur.  For  the  case  u+  —  0.1,  the  special  value  is  approximately 
=  0.332051.  A  detailed  discussion  of  the  phase  portrait  for  this  example  is  given  in 
[BMS99].  Shown  in  Figure  4  is  a  Poincare  section  (at  fixed  u)  from  the  phase  portrait 
of  (2.30).  The  stable  manifold  WS(B)  is  shown  as  a  dashed  line  while  the  unstable 
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Figure  4:  Section  of  the  invariant  manifolds  WU(M )  (solid  line),  WS(B)  (dashed  line)  and 
WU(T )  (plus)  with  the  Poincare  plane  P  —  {(it,  u' ,  it");  it  =  (2 it_  +  tt+)/3).  Intersections 
of  the  solid  and  dashed  lines  mark  points  where  trajectories  connecting  M  to  B  hit  P. 
The  symbols  around  these  intersections  indicate  whether  the  corresponding  profile  is  a 
stable  (circles)  or  unstable  (boxes)  solution  of  the  PDE,  for  the  first  two  traveling  waves 
u.  Arrows  indicate  the  direction  of  ((<p^),  (‘ft)')  and  ((^3  );  (‘fs)')  (the  projection  of  the 
derivative  vectors  on  the  Poincare  plane). 


manifold  WU(M )  is  shown  as  a  solid  line.  Note  that  this  unstable  manifold  appears  as  a 
spiral  whose  center  corresponds  to  the  undercompressive  connection.  Each  intersection 
of  these  two  curves  denotes  a  point  along  a  heteroclinic  orbit  of  (2.30).  The  intersection 
denoted  by  a  circle  corresponds  to  a  one-dimensionally  stable  compressive  wave  while  the 
box  corresponds  to  a  one  dimensionally  unstable  compressive  wave.  The  arrows  on  the 
figure  denote  the  orientation  of  a  choice  of  basis  vectors  (<p{~),  (‘ft)')  and  ((^3  ),  (^3)')- 
Note  that  the  orientation  of  these  vectors  switches  at  each  successive  crossing  of  the 
spiral.  The  vector  (ux,  u'x,  u")  is  into  the  Poincare  plane.  Following  Corollary  2.8,  the 
corresponding  stability  indices  Tj  alternate  in  sign.  This  is  consistent  with  numerical 
observations  that  the  successive  waves  alternate  stability.  (Recall  that  u±  are  fixed.  The 

,  /  ux  ipt  \  1  ,  /  ux  \  .  . 

remaining  terms  det  7/  and  det  _#  fix  a  common  orientation  on  the 

V  UX  <P2  )  V  ^3  Ux  J 

stable/unstable  manifolds  at  u_/u+). 

Remark.  As  often  happens,  we  obtain  a  rigorous  instability  result  by  this  tech¬ 
nique,  but  it  is  only  suggestive  of  stability.  To  obtain  a  complete  stability  result  would 
require  establishing  nonexistence  of  pure  imaginary  eigenvalues  (other  than  A  =  0),  a 
separate  (and  quite  interesting)  unresolved  issue. 
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2.6  Evaluation  at  A  =  0  in  the  undercompressive  case. 

The  undercompressive  case,  sgn  a_  =  sgn  a+,  can  be  treated  similarly.  For  example, 
Proposition  2.6  goes  through  unchanged.  However,  the  behavior  at  A  =  0  is  significantly 
different,  as  might  be  expected  from  the  fact  that  long- wave  behavior  is  dominated  by 
the  convection  rates  a±  in  the  far  held.  In  the  present,  one-dimensional  context,  this 
leads  to  a  slightly  modified  stability  index.  In  the  context  of  multi-dimensional  stability, 
as  we  shall  see  later,  the  distinction  is  still  more  critical. 

Without  loss  of  generality  (and  to  be  consistent  with  the  thin  him  context  [BMS99]), 
assume 


a-  =  /'(«-),  a+  =  f\u+)  <  0. 


(2.37) 


Note  that  this  implies  that  there  exists  an  intermediate  equilibrium  5  um  between 
and  u. |_.  The  corresponding  system  of  equations  (2.30)  describing  the  phase  portrait 
for  the  traveling  wave  has  (at  least)  three  equilibria  B  —  (u+,  0,0),  M  —  (um,  0,0), 
T  —  (it_,0,  0).  The  undercompressive  wave  corresponds  to  a  connection  from  T  (which 
has  a  one  dimensional  unstable  manifold)  to  B  (with  a  two  dimensional  stable  manifold). 

Appealing  again  to  Section  2.2,  we  hnd  that,  as  in  the  Lax  case,  <S+  consists  of  all 
solutions  of  (2.5)  decaying  exponentially  at  x  —  +oc.  However,  recalling  the  discussion 
before  Corollary  2.2,  when  A  =  0,  U~  now  consists  of  solutions  that  are  only  bounded  as 
x  — >•  — oo.  Without  loss  of  generality,  when  A  =  0,  choose  again  <Pi(x)  —  ipf(x)  —  ux(x), 
and  take  tpf,  cpf  to  be  independent  solutions  of  (2.5)  in  <S+,  IA~  respectively,  where  we 
can  choose  the  normalization  ip 3  (— oo)  =  1. 


Proposition  2.9  D( 0)  =  0,  while 


<P2'  (z  ~ z+) 

sgn  D'(0)  =  -  sgn  (a_/c(0))det  (  u'x  ip%'  (z~  -  z+)'  |  |I=0 

Ux  Vz  (■ z  ~Z+) 


u.r 


=  sSn  a-  JZc([/c)(u  -  u-)  det  (  J 


dx, 


(2.38) 


where  satisfies  the  variational  equation 


cz"  —  hz!  —  az  —  (u  —  U-) , 


(2.39) 


with 


z+(+oo)  =  0,  z  (— 00)  =  (u+  —  w_)/a_. 

5Note  that  we  are  using  slightly  different  notation  from  [BMS99]  in  which  u_ 
nearest  equilibria. 


(2.40) 

and  u+  always  denote 
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Proof.  As  in  the  proof  of  Proposition  2.7, 


u., 


Pi 


P3  ux 


\ 


D(  0)  =  det 


=  0, 


(2.41) 


and 


D'(  0)  =  det 


<  {piT  (p3)'"  <  ) 

<  ux  p+  (f3  ^ 

\  ux  pj  p3  \P  4\  -<PiJ  ) 


(2.42) 


Likewise,  the  decaying  solutions  ux  and  p \  satisfy  equation  (2.28).  However,  p3 
is  asymptotically  constant,  p3  (— 00)  =  1,  hence  satisbes  the  inhomogeneous  equation 


op"'  —  bp1  —  cap  +  a_. 


(2.43) 


Using  the  fact  that  z  :  =  <p4  —  ipf  as  before  satisbes  (2.35),  we  can  use  (2.28),  (2.35), 
(2.43)  to  reduce  the  fourth  row  of  (2.42),  obtaining 


D'(  0)  —  c  1  det 


/  ux  <p$  p3,  z  \ 

K  pi  p3  % 

— //  -\-n  ~if 

Ux  Pi  P3  z 

\  0  0  a_  u_  —  u+  ) 


(2.44) 


Using  the  third  column  to  eliminate  (u+  -  m_),  we  obtain 
D'(  0)  =  c_1  det 


f 

^ X 

pi 

p3/ 

z 

-  z+ 

K 

pi 

pi 

(z~ 

-  z+y 

n 

Ux 

pi" 

n 

P3 

(z- 

-  z+y 

V 

0 

0 

a_ 

0 

(2.45) 


where  2  —  z+  —  z  — 


U—  —U-\- 

a 


p3  are  debned  by 


Pl\ 


U_  —  U. |_ 


a_ 


P3 


(2.46) 


Evaluation  of  (2.45)  at  x  —  0  gives  the  brst  equality  in  (2.38).  Integration  of  (2.34), 
together  with  (2.43),  verify  that  z±  both  satisfy  (2.39),  with  the  claimed  boundary 


conditions  at  ±00. 

Now  debne 

/  ux 

pi 

z  \ 

4/±  ;=  det  (  u'x 

pi 

z+'  . 

(2.47) 

\< 

pi" 

z+"  ) 
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The  reduction  to  a  Melnikov  integral  in  the  second  equality  in  (2.38)  involves  evalua¬ 
tion  of  —  4/+.  Briefly,  this  follows  (see,  e.g.,  [GZ],  proof  of  Lemma  3.4,  or  [GH,  Sch]) 
from  the  (inhomogeneous)  Abel’s  formula, 


ip1  —  TrM  ip  +  det 
f  u*  Vi 

—  det 


/  Ux  (p 2  0  ^ 

:  0 

vt"  («-«-)/ 

°  \ 

0 

V  K  V2"  (u  -  «-)  / 


(2.48) 


and  Duhamel’s  principle/ variation  of  constants,  to  evaluate  4,±  at  x  —  0.  Here,  B  = 

/  0  1  0\ 

1  denotes  the  coefficient  matrix  for  (2.28)  written  as  a  first  order 


,x) 


—a/c  b/c  0 


\  / 

system.  ■ 

Combining  Corollary  2.8  and  Proposition  2.9,  and  recalling  that,  as  x  — >•  —  oo, 
(03 ,03  )  -»  (1,0)  and  sgn  uxx  ~  sgn  ^ ux  =  sgn  ux,  we  have 


Corollary  2.10  The  stability  index  T  :=  sgn  D'(0)D(+oc)  is  given  by 


r  = 


sgn  a_ 


)  det 


ux  pi; 

K  vt 


dx  ux(—oo)  det 


ux  (p% 

-i  +' 

ux  Vi 


£  =  +  00 
(2.49) 


Note  that  the  functions  z±  in  (2.38)  are  uniquely  specified  by  (2.39)-(2.40),  modulo 
spanllL:,  <^}.  and  z+  have  an  important  interpretation  as  the  variations,  in  the 
traveling  wave  ODE, 


c{u)u"  —  b{u)u  —  ( f(u )  —  f(u-)  —  s(u  —  ii_)),  (2.50) 


of  the  unstable/stable  manifolds  of  u+  and  w_(s)  as  the  shock  speed  s  is  varied.  Likewise, 
(  u*  vi  (z~  -  Z+)  \ 


:(0) 


det 


|cc=0 


(1  /c)(u  —  u_)  det 


\  ux  V2  (z  ~  z+)"  ) 


Ux  (p  2 

K  vt' 


(2.51) 


can  be  recognized  as  |j|s=0,  where  d(s)  is  the  Melnikov  separation  function,  measur¬ 
ing  the  distance  (in  3-dimensional  phase  space)  between  the  one  dimensional  unstable 
manifold  T  and  the  two  dimensional  stable  manifold  of  B  along  a  fixed  transverse  sec¬ 
tion.  Alternatively,  it  represents  the  orientation  with  which  <S+  and  U~  intersect  in 
4-dimensional  phase  space  of  eqn.  (2.50)  augmented  with 


s  =  0. 


(2.52) 
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The  stability  index  T  is  a  coordinate-  and  normalization-independent  measure  of  the 
orientation. 

Calculation  of  T. 

We  now  turn  to  the  calculation  of  the  stability  index  for  the  undercompressive  case, 
which  hinges  on  the  calculation  of  the  sign  of  the  Melnikov  integral 

dd/ds—  f  (l/c)(u  —  U-)  det  f  Jf  ^  dx. 

J  -  oo  \Ux  ¥ 2  / 

This  can  be  evaluated  conveniently  in  terms  of  the  left  zero  eigenfunction  7r°.  We  post¬ 
pone  a  discussion  of  the  associated  spectral  theory  to  the  next  section.  It  transpires  that 
the  Melnikov  integral  arising  in  calculation  of  the  stability  index  for  undercompressive 
shocks  can  be  expressed  alternatively  as 

{u  —  it*,  kn®},  —  (ux,  kn0)  —  k,  (2.53) 

where  (•,  •)  denotes  L 2  inner  product,  u*  is  the  state  arising  in  the  original  formula  (in  the 
present  setting,  u*  —  w_;  see  [GZ]  for  more  general  situations),  and  k  is  an  appropriately 
chosen  constant.  The  key  simplification  of  the  Melnikov  integral  comes  from  the  formula 

kTT°x  =  (1/c)  det  f  ^  ,  (2.54) 

which  we  derive  below.  Comparing  behaviors  at  +oo,  we  find  that  the  sign  of  the 
Melnikov  integral  is 

sgn  (1/c)  det  ^|/^)  tt° |i=+00, 

and  thus,  referring  back  to  (2.50),  we  obtain  the  simple  expression: 

Proposition  2.11  For  an  undercompressive  wave,  the  stability  index  T  can  be  computed 
as 


T  =  sgn  a_7r^(+oo)Mx,(— oo)  (2.55) 

where  7ru  denotes  the  left  zero  eigenfunction. 

For  the  specific  traveling  wave  computed  numerically  and  shown  in  Figure  5  and 
for  the  corresponding  7r°  shown  in  Figure  7,  we  observe  that  T  =  +1,  consistent  with 
stability  of  the  undercompressive  wave.  To  finish  the  proof  of  Proposition  2.11  we  show 
the  equality  (2.54)  then  use  integration  by  parts  and  the  fact  that  7r°(+oo)  =  0.  As 
we  discuss  in  more  detail  in  the  next  section,  7ru  is  uniquely  specified  up  to  a  constant 
factor,  is  bounded,  and  satisfies  the  adjoint  eigenvalue  equation 

—  (cz1)111  +  (bz1)'  +  az  —  0. 
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It  follows  that  (ei,e2,e3)  must  satisfy  the  adjoint  ODE  to  the  linearized  traveling 
wave  equation  written  as  first  order  system  W'  —  BfW,  where  13  is  as  in  (2.6), 


from  which  we  readily  find  that  e3  satisfies 

e3  -  (be3/c)'  -  ( a/c)e3  =  0. 


Thus, 

(l/c)detf“f  %\={\!c)ez 

\  Ux  r 2  J 

satisfies  (2.56)  as  claimed. 

Remarks. 

The  formula  (2.55)  in  Proposition  2.11  is  new.  The  improved  formula  (2.56)  in 
Proposition  2.11  is  not  restricted  to  the  special  setting  of  thin  film  models.  In  fact,  the 
ideas  used  in  its  derivation,  in  particular  (2.53)  and  (2.54),  easily  generalize  to  systems 
of  conservation  laws  of  arbitrary  order  to  yield  an  analogous  result.  These  formulae  are 
the  consequence  of  a  general  but  somewhat  nonstandard  duality  principle.  6 

This  observation  therefore  represents  a  useful  refinement  in  the  general  theory.  For 
comparison,  the  formula  obtained  in  [GZ]  for  second  order  diffusive  systems  is  analogous 

°In  particular,  note  that  the  fact  that  e^/c  satisfies  the  equation  for  z'  can  be  seen  from  the  more 
general  relation 

(z' ,  z" ,  z"  )S(w,  w' ,  w"  Y  =  constant 

for  solutions  w  of  the  linearized  traveling  wave  ODE  (in  particular,  decaying  solutions  of  the  adjoint 
ODE  at  A  =  0)  and  solutions  z  of  the  adjoint  eigenvalue  ODE  at  A  =  0,  where 

f-b  +  C"  -J  <\ 

S  :=  2c'  -c  0  .  (2.57) 

\  c  0  0/ 

For,  by  duality,  this  implies  that  (z1 ,  z" ,  z"')S  =:  (ei,e2,es),  solves  the  adjoint  of  the  traveling  wave 
ODE  written  as  a  first  order  system,  and  by  inspection  e-$  =  cz' . 

This  relation,  in  turn,  is  a  special  case  of 

[z,z,z  ,z  )o(w,w  ,w  ,w  )  =  constant 

for  solutions  w,  z  of  the  eigenvalue  and  (resp.)  adjoint  eigenvalue  ODE  at  A  =  0,  where 

(a  b  +  c"  —c'  c\ 

-b-c"  2c'  c  0  ,  (2.58) 

-c'  -c  0  0/ 

under  assumption  cw'"  =  bw'  —  aw.  The  duality  relation  (2.57)  was  pointed  out  in  [LZ.2]  in  the  second- 
order  diffusive  (system)  case;  the  underlying  relation  (2.58)  was  pointed  out  in  [ZH],  Lemma  4.4,  and 
the  extension  to  arbitrary  higher  order  operators  discussed  in  [ZH],  proof  of  Theorem  6.3. 
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to  Corollary  2.10;  in  the  third  order  scalar  case  treated  in  [D]  the  formula  obtained  are 
explicitly  evaluable,  and  the  issue  does  not  arise. 

As  pointed  out  in  [SSch]  in  the  diffusive  case,  nonvanishing  of  the  Melnikov  integral 
is  related  to  well-posedness  of  Riemann  solutions  near  the  undercompressive  shock  (w_ , 
u. |_,  s).  There  is,  however,  no  analog  of  (2.31)  in  the  diffusive  case,  for  which  scalar  Lax 
shocks  are  always  stable  [S.l,  H.l,  H.2,  H.3,  JGK]. 

In  both  Lax  and  undercompressive  cases,  the  index  I  can  be  evaluated  numerically. 
However,  it  does  not  give  a  conclusive  result  of  stability,  yielding  only  parity  of  the 
number  of  unstable  eigenvalues  of  L. 

A  complete  numerical  evaluation  of  stability  can  be  performed  efficiently  using  an 
algorithm  of  Brin  [Br],  in  which  D(-)  itself  is  evaluated  numerically,  and  the  winding 
number  computed  around  {.Re A  >  0).  Moreover,  the  top  eigenvalue  can  be  computed 
numerically  using  the  power  method  discussed  in  the  following  section.  Furthermore,  as 
was  done  in  [BMS99],  nonlinear  stability  can  be  verified  numerically  by  perturbing  the 
stationary  wave  and  noting  convergence  or  divergence  of  solutions  of  the  PDE. 

3  Multi- dimensional  stability  and  the  long-wave 
paradox. 

This  section  addresses  multi-dimensional  stability  and  the  long-wave  paradox  for  stabil¬ 
ity  of  undercompressive  shocks.  Here  we  are  concerned  with  how  the  spectrum  of  the 
linearized  operator  depends  on  the  transverse  wave  number.  For  this  problem  we  present 
theoretical  justification  for  formal  and  numerical  treatments  of  the  spectral  problem. 

A  first  order  perturbation  analysis  of  the  spectrum  can  be  readily  carried  out,  both  for 
systems  and  for  scalar  equations,  by  Evans  function  calculations  very  similar  to  those  of 
the  previous  section  [ZS] .  What  we  require  here,  however,  is  a  second  order  expansion, 
or  ‘diffusive  correction’  in  the  language  of  [ZS],  since  scalar  shock  fronts  are  always 
neutrally  stable  to  first  order  (see  [ZS],  or  calculations  just  below).  In  the  present,  scalar 
setting,  this  is  much  more  convenient  to  carry  out  via  standard,  Fredholm  solvability 
calculations,  justified  by  passing  to  an  appropriate  weighted  norm,  than  via  the  Evans 
function  framework,  and  we  shall  therefore  take  this  different  point  of  view  in  what 
follows. 

3.1  Preliminaries 

Consider  now  the  multi-dimensional  version 

ut  +  f{u)x  =  V  •  (6(m)V-u)  -  V  •  (c(w)VArt)  (3.1) 

of  equation  (2.1).  Linearizing  about  u(-),  and  taking  the  Fourier  transform  in  transverse 
directions,  we  obtain 

vt  =  Lkv  :=  ~(cvxxx)x  +  (bvx)x  -  ( av)x  -  k2(bv  -  ( cvx)x  -  cvxx)  -  k^cv 

(3.2) 
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where 


c  =  c(u),  b  =  b(u),  a  —  f(u)  —  b'(u)ux  +  c'(u)uxxx,  (3.3) 

and  k  is  the  wave  number,  v  —  v(x,  k,t). 

Remark.  Equation  (3.2)  is  for  two  space  dimensions  as  in  the  application  to  thin 
films.  In  the  case  of  higher  dimensions,  rotational  invariance  yields  the  same  transformed 
equation  in  which  k  becomes  the  modulus  of  the  wavenumber. 

Associated  with  (3.2),  we  have  the  family 

(L/,  —  X  I)iu  —  0  (3.4) 

of  eigenvalue  equations,  where  Lq  is  exactly  the  linearized  operator  about  u(x)  for  the 
one-dimensional  problem.  In  what  follows,  we  will  assume  L ^  is  an  analytic  function  of 
k2,  rather  than  just  k,  since  in  equation  (3.2),  L ^  is  a  quadratic  polynomial  of  k2. 

The  question,  assuming  stability  of  the  one-dimensional  operator,  is  whether  the 
spectrum  of  L &  remains  stable  as  k  is  varied.  Of  particular  interest  (see  e.g.  [CE]) 
is  the  variation  for  small  k  of  the  top  eigenvalue  X(k),  where  A(0)  =  0  corresponds  to 
translation  of  the  wave.  For  it  is  this  mode,  to  leading  order,  that  governs  the  propagation 
of  disturbances  along  the  front  (see  discussion  [G,  GM,  HZ. 2,  K.l,  K.2,  K.3,  KS,  ZH,  ZS] ) . 
Stability  of  X(k)  for  small  k  (i.e.  Re X  <  0)  is  called  long-wave  stability.  This  is  the 
problem  of  most  interest  since  the  fourth  order  diffusion  necessarily  causes  the  dominant 
eigenvalue  to  decay  as  —k4c  for  large  values  of  \k\. 

3.2  Spectral  Theory. 

A  rigorous  discussion  of  long- wave  stability,  in  the  context  of  conservation  laws,  requires 
additional  spectral  perturbation  theory,  beyond  the  usual  Banach  space  theory  of,  e.g.  [K, 
Y].  For,  the  eigenvalue  A(0)  =  0  of  main  interest  is  embedded  in  the  essential  spectrum 
of  Lq,  [He,  S.l,  ZH]  hence  this  standard  theory  does  not  apply.  A  priori  there  is  no 
reason  to  expect  an  analytic  development  of  A (k)  at  such  a  point;  indeed,  for  systems  of 
conservation  laws,  it  is  a  fundamental  fact  that  A(-)  is  in  general  not  analytic  at  k  —  0. 
We  refer  the  reader  to  [ZS]  for  further  discussion  of  this  general  situation. 

However,  in  the  present,  scalar  context,  a  much  simpler  treatment  is  possible  by  the 
weighted  norm  method  of  Sattinger  [S.2] .  Introducing  the  norm 

ll/l|n:=ll/S2|li>,  (3.5) 


with 


n(x)  :=  e-ef°a{y)dv,  1  >  0  >  0,  (3.6) 

has  the  effect  of  shifting  the  essential  spectrum  of  Lq  to  the  left,  into  the  strictly  stable 
complex  half-plane,  as  is  readily  checked  by  the  methods  of  [He,  S.2]  (see  also  [Z.2]  for 
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further  details).  Thus,  in  this  norm,  A(0)  =  0  becomes  an  isolated  eigenvalue,  and  we 
may  conclude  from  standard  spectral  theory  the  existence  of  analytical  developments 

A  (k)  =  A1  A;2  +  •  •  •  ,  (3.7) 

and 

<p(k)  —  <p°  +  tp1k2  H - ,  (3.8) 

of  A  and  the  associated  right  eigenfunction  ip  £  Lq. 

Likewise,  there  exists  an  analytic  L q  left  eigenfunction  with  respect  to  the  Lq  inner 
product,  whence  we  can  conclude  existence  of  an  analytic  eigenfunction 

7 v(k)  —  7T°  +  7 r1^2  +  •  •  •  ,  (3-9) 

with  respect  to  the  standard  inner  product,  where  f2_27r  G  or  equivalently 

vr  ELl-i.  (3.10) 

Moreover,  A / <p  are  the  unique  Lq  eigenvalue/function  pair  for  in  the  vicinity  of  A  =  0, 
and  A/7 r  the  unique  1  pair. 

An  interesting  feature  of  these  developments  as  compared  to  the  usual  (unweighted) 
type  is  that,  depending  on  the  signs  of  a±,  p(k)  or  n(k)  may  be  nondecaying,  in  fact 
exponentially  growing  at  ±00.  Nonetheless,  as  in  the  classical  case,  we  have  the  following 
result,  in  which  (•,  •)  denotes  the  standard  L2(M)  inner  product: 

Proposition  3.1  Suppose  that  (consistent  with  one- dimensional  stability)  X  —  0  is  a 
simple  eigenvalue  of  Lq  (i.e.,  when  k  —  0),  and  that  all  other  eigenvalues  have  negative 
real  part.  Let  A (k)  as  above  denote  the  top  eigenvalue  of  L Then,  for  compactly 
supported  initial  data  vq,  and  x  restricted  to  any  bounded  domain,  the  solution  v(x,t)  — 
eLktVo(x)  of  vt  —  LkV  satisfies 

v(x,t)  =  ex^tT{k)v 0(1  +  0(e-^)),  p  >  0,  (3.11) 

where  V{k)  denotes  the  generalized  spectral  projection  operator 

'P(k)g  :=  p(k){n(k),g)  (3.12) 


associated  with  X(k). 

This  can  be  seen  from  the  corresponding  classical  result  in  L^,  together  with  the 
observation  that  is  bounded  above  and  below  on  compact  domains.  For  a  proof  in  the 
general,  systems  setting,  see  [ZH],  sections  5-6  and  8  (note:  this  result  is  independent  of 
the  regularity  of  A(-).) 

Proposition  3.1  states  that  the  nonstandard  eigenfunctions  ip,  7r  still  govern  asymp¬ 
totic  behavior  on  bounded  domains,  analogous  to  “resonant  poles”  in  scattering  theory 
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[LP] .  A  useful  corollary  of  this  fact  is  that  the  power  method  can  be  used  to  find  the 
top  eigenvalue/eigenfunction  of  Lk,  by  solving  vt  —  L^v  starting  with  any  compactly  sup¬ 
ported  initial  data.  This  validates  the  numerical  stability  analysis  of  [BB97,  BMFC98]. 
Likewise,  the  classical  formula 

A1  =  {7T°,  LV>  +  (7T°,  (L°  —  A0)^1)  =  {7T°,  LV)  (3.13) 

(Lk  L°  +  kL 1  +  •  •  • )  based  on  solvability  conditions/  Fredholm  alternative  can  be 
seen  to  apply,  as  a  consequence  of  the  equivalent  formula  in  the  f2  norm.  That  all  L2 
inner  products  are  well-defined  follows  from  consideration  of  the  decay/growth  of  ip,  7r 
as  x  — y  ±oo,  which  reveals  that  they  and  their  derivatives  lie  in  Lq,  respectively. 

(By  comparison,  only  a  weak  version  of  the  Fredholm  alternative  survives  for  systems, 
see  [ZH]). 

Remark.  The  introduction  of  the  weighted  norm  ||  •  ||n  is  equivalent  to  the  change  of 
variables  z  —  Qw,  from  which  observation  we  may  deduce  that  (since  such  a  change 
of  variables  introduces  only  a  nonzero  constant  multiplier)  that  zeroes  of  the  Evans 
function  correspond  to  eigenvalues  with  respect  to  f2  norm.  Moreover,  this  yields  the 
useful  characterization  of  their  generalized  eigenfunctions  as  the  intersection  between  the 
stable/unstable  manifolds  S+  /U~  defined  as  in  Section  2. 

To  investigate  long-wave  stability,  we  expand  the  operator  Lk  in  powers  of  k2,  Lk  — 
L°  +  ArL1  +  A;4L2,  with 

£l  =  -b  +  ^  (CK)  +  ci~I-  £2  =  -c-  (3.14) 

Using  the  fact  that  ip°  —  ux,  to  first  order  we  obtain 

A1^,  =  Ll7ux  +  L^cp1. 


We  multiply  this  equation  with  the  left  eigenfunction  of  Lq,  i.  e.,  the  solution  of 


L*y  =  o, 


(3.15) 


where 


is  the  formal  adjoint,  and  then  integrate  from  —  oc  to  +oo,  to  get 

/OO  /»oo  /»oo 

n°ux  —  /  Tr0^1^.  +  /  7 r°L°(/)1. 

■oo  J  —  OO  J  —  OO 


T*  — 

L  o  — 


d 3 


dx3 


d  \  d 
dx  J  dx 


Integration  by  parts  of  the  last  term  yields 


7 r°Ly 


klL*0  7T°  =  0, 


(3.16) 


(3.17) 
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since  the  boundary  terms  vanish. 

Hence, 

/OO  poo 

ty°ux  =  /  nQl}ux, 

•oo  J  —  OO 

as  in  (3.13)  if  7T°  is  normalized  by  <  7 r°,ux  >—  1. 

Hence,  using  the  explicit  expression  for  L1  in  (3.14),  we  obtain,  after  integration  by 
parts 

/oo  poo  poo 

7T°ux  =  7T°ct2a;a;|“oo  -  /  n°(bux  -  cuxxx )  -  /  TT^,CUXX . 

•oo  J  — OO  J  — OO 

Since  the  boundary  terms  vanish,  the  traveling  wave  ODE  implies 


/OO  /»oo  /»oo 

=  -  /  7T°(/(iZ)  -  /(«_))  -  /  7T°ciZ; 

•oo  J  —  oo  J —oo 


(3.18) 


To  proceed  further,  we  have  to  determine  the  left  eigenfunction,  by  finding  a  solution 
of  (3.15)  with  the  appropriate  behavior  at  ±oo. 

The  discussion  of  the  boundary  conditions  for  7T°  depends  on  the  growth/decay  rate 
of  the  traveling  wave  and  the  first  order  correction  of  the  left  eigenfunction,  and  leads 
to  different  results  for  the  compressive  and  the  undercompressive  case. 

In  the  Lax  or  compressive  case,  (a_  >  0  >  a+),  we  have  that  the  weight  Q  of 
(3.6)  grows  exponentially.  Thus,  functions  in  Lq  decay  exponentially  at  Too/  — oo,  while 
functions  in  T^-i  are  allowed  to  it  grow  exponentially.  Thus,  the  constant  solution 


L*0w  =  ~{cwx)xxx  +  (bwx)x  +  awx  =  0  (3.19) 

is  the  unique  left  eigenfunction,  whereas  ux  is  the  exponentially  decaying  right 

eigenfunction  in  L fl.  (Note:  (7r0,!4,)  =  1).  This  validates  the  formal  calculation  in 
[THSJ89,  BB97]  via  (3.13),  which  in  this  case  is  equivalent  to  simple  integration,  yielding 

A1  = 

In  the  undercompressive  case,  (a_,a+  <  0),  on  the  other  hand,  Q  blows  up  as 
x  — >  Too  but  decays  as  x  — >  — oo.  Thus,  L q  functions  may  grow  exponentially  while 
functions  must  decay  exponentially  at  — oo.  The  latter  fact  implies  that  7r°  in  this 
case  is  not  a  constant  function,  hence  the  formal  calculation  of  [THSJ89,  BB97]  no  longer 
corresponds  to  (3.13). 

Moreover,  numerical  computations  by  the  power  method  reveal  that  ip1  indeed  grows 
exponentially  as  x  — >  — oo.  Thus,  the  calculation  of  the  Lax  case,  by  integration  against 
a  constant  is  not  valid  to  test  stability  in  the  undercompressive  case,  the  resulting 
integral  being  unbounded,  hence  the  conclusion  of  instability  by  this  test  is  false.  (Of 


U. |_  —  U— 


dx. 


(3.20) 
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course,  this  fallacy  would  likewise  be  detected  at  the  level  of  solving  for  ip1.)  This  gives 
a  simple  explanation  of  the  aforementioned  long- wave  paradox. 

We  further  examine  the  growth  rate  of  p> 1  at  —  oo  analytically:  The  slow  growth 
rate  p-~(k,  A)  perturbing  from  the  one-dimensional  root  /i_(0,  0)  =  0  of  the  characteristic 
equation  for  Lk±, 


—c±/ji  +  b±p2  —  /j,a±  —  c±k 4  —  k2(b±  —  2c±p 2)  =  A 

(analogous  to  the  one-dimensional  version  (2.7)),  at  k  —  0,  A  =  0,  gives 

(A  +  b±k 2  +  c±k 4) 

p±(k,  A)  = - . 

a± 


(3.21) 


(3.22) 


Now,  let  us  consider  the  case  6  =  0  corresponding  to  the  stability  computations  of 
the  undercompressive  wave  in  [BMFC98].  Here,  /i_  =  —  A/a_  +  H.O.T..  Thus,  in  case 
of  stability,  A1  <  0,  we  have 


P-(k,  \(k))  ~  p-(k,  A1  A;2)  ~ 


A1*;2 

- <0, 

a_ 


(3.23) 


and  so  p(k)  generically  exhibits  exponential  growth  as  x  — >  — oo  for  any  small  k  >  0, 
the  exception  arising  when  <S+  intersects  U~  precisely  along  the  unique  decaying  solution 
corresponding  to  the  remaining  (positive)  root.  Indeed,  ip(k)  (generically)  decays  at  ±oc 
in  case  of  instability! 


Remark.  The  growth  of  p  at  —  oo  raises  also  an  apparently  subtle  issue  of  competition 
between  exponential  temporal  decay  eReX(k)1  and  exponential  spatial  growth,  and  indeed 
this  issue  is  not  resolved  by  considerations  in  the  f2  norm  (where  spatial  growth  is  ignored 
at  — oo).  However,  heuristically,  note  that  at  the  order  of  Ai,  the  dominant  linear  behavior 
in  the  far  held  at  —  oo  is 

ek2\1te~\1k2x/a-  =  exp(lMl(x  _  a_fy 

which  is  precisely  a  translating  wave  with  speed  a_  and  shape  corresponding  to  the 
eigenfunction  in  the  far  held.  Thus  the  exponential  growth  of  the  eigenfunction  and 
exponential  decay  of  the  eigenvalue  combine  to  produce  a  translation  in  the  far  held. 
This  is  consistent  with  the  fact  that  for  the  undercompressive  wave,  information  passes 
through  the  waves  and  goes  oh  to  —  oo  with  speed  a_.  Preliminary  numerical  compu¬ 
tations  [M99]  verify  that,  in  this  regime,  temporal  decay  indeed  wins  out  over  spatial 
growth. 

More  detailed  computations  of  [HoZ.3,  Z.l],  in  the  second  order  diffusive  case  carry 
over  in  straightforward  fashion  to  the  case  of  higher-order  diffusion  (see  [HZ. 2]  for  a 
one-dimensional  version)  to  reveal  that 

( Vk )  ReX(k )  <  9k2,  9  <  0,  for  the  top  eigenvalue  X(k)  of 
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is  a  sufficient  condition  for  linearized  stability,  independent  of  spatial  growth  of  eigen¬ 
function  if,  or  in  our  notation, 

A1  <  0.  (3.24) 


3.3  Alternative  solvability  condition. 

As  discussed  above,  in  the  undercompressive  case,  7T°  =  constant  is  no  longer  a  left 
eigenfunction  for  ux  (as  in  the  Lax  case).  Indeed,  this  is  the  key  difference  between  Lax 
and  undercompressive  waves  pointed  out  in  [LZ.2] :  the  shock  shift  (represented  at  the 
linearized  level  by  the  component  of  instantaneous  translation  ux )  is  not  determined  by 
mass,  (i.e.  projection  (1,-)),  but  by  a  different  time-invariant  of  the  solution  (7T° (•) ,  •), 
where  7T°  is  an  nonconstant  bounded  solution  of 

L*ir°  =  0,  (3.25) 


7r°(a:)  — >  1 
7r°(a:)  — >■  0 


as  x  — >•  +oc 
as  x  — >•  —  oo. 


(3.26) 


The  boundary  conditions  can  be  understood  heuristically  from  the  observation  that 
mass  near  Too  is  swept  inward  by  convection  a+  <  0,  to  interact  with  the  shock,  while 
mass  near  — oo  is  swept  away  by  a_  <  0,  never  affecting  the  shock  (see  [ZH],  §10  for  a 
rigorous  discussion).  The  condition  at  Too  may  be  deduced  rigorously  by  consideration 
of  growth/decay  rates  at  +oo  of  solutions  of  the  adjoint  eigenvalue  equation,  which 
reveals  that  solutions  in  must  in  fact  be  bounded  (note:  growing  solutions  grow  at 
too  great  an  exponential  rate).  The  condition  at  +oc  is  clear,  by  exponential  decay  of 
functions  in  . 


3.4  Numerical  Results 

In  this  subsection,  we  describe  numerical  results  for  the  thin  film  equation 

ut  +  ( u 2  —  u3)  x  —  — V  •  (ti3VAu)  ,  (3.27) 

which  is  equation  (1.5)  with  b  —  0,  c  —  u3,  f  —  u2  —  u3.  We  focus  exclusively  on 
multidimensional  stability  of  one-dimensional  traveling  waves.  A  preliminary  stability 
analysis  in  [BMFC98]  of  compressive  waves  revealed  that  they  were  unstable  against 
2-dimensional  perturbations  for  a  range  of  wavenumbers  k  >  0.  On  the  other  hand, 
the  undercompressive  wave  was  stable  for  all  wavenumbers.  These  results  agree  with 
experimental  observations  of  the  liquid  front. 

In  this  paper,  we  calculate  the  critical  growth  rate  A (k)  for  a  variety  of  compres¬ 
sive  and  undercompressive  traveling  waves,  and  compare  the  results  with  the  long-wave 
asymptotic  result  A (k)  ~  A1^2  as  k  — >  0+,  in  which  the  coefficient  is  calculated  numer¬ 
ically  from  the  formulae  above.  There  are  several  interesting  features  of  the  behavior  of 
A (k)  for  moderate  k,  linked  to  the  presence  of  a  countable  family  of  compressive  waves 
for  the  wave  speed  at  which  there  is  also  an  undercompressive  wave. 
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3.4.1  Computation  of  the  largest  eigenvalue  using  the  power  method 

For  each  of  the  traveling  wave  solutions,  the  values  of  the  base  profiles  u  on  the  finite 
difference  grid  are  obtained  by  solving  the  traveling  wave  ODE,  or  by  computing  long¬ 
time  solutions  of  the  one-dimensional  PDE.  Then  the  dominant  eigenvalue  X(k)  and  the 
corresponding  eigenfunction  are  obtained  for  various  wavenumbers  k  by  solving  (3.2) 
numerically.  As  in  [BB97,  BMFC98],  we  used  a  finite  difference  scheme,  with  an  implicit 
Euler  time  step  to  calculate  v,  then  extracted  the  eigenvalue  from  the  time  evolution 
of  v  after  the  exponential  decay  or  growth  rate  in  time  had  appeared.  In  addition, 
eigenvalues/functions  were  also  confirmed  by  inverse  vector  iteration,  a  particularly  use¬ 
ful  method  for  computing  subdominant  or  pairs  of  complex  conjugate  eigenvalues  as  in 
section  3.4.3. 

We  set  parameters  for  the  travelling  waves  by  choosing  u+  —  0.1,  U-  —  0.332051, 
and  shift  to  a  reference  frame  (by  subtracting  a  linear  term  from  the  flux)  so  that  the 
wave  is  stationary.  For  this  choice  of  parameters,  we  find  an  infinite  number  of  traveling 
wave  solutions,  as  observed  in  [BMS99],  but,  as  discussed  in  Section  2.5,  only  every 
other  one  is  stable  as  a  solution  of  the  one- dimensional  version  of  (3.27).  For  the  same 
parameters,  there  are  two  additional  travelling  waves:  an  undercompressive  wave,  from 
Ut  =  1  —  it-  —  u+  —  0.567949  to  u+,  and  a  trailing  compressive  wave  from  to  Ut- 
(The  latter  wave  is  termed  trailing  because  it  corresponds  to  compressive  waves  that  trail 
the  undercompressive  wave  in  numerical  simulations  [BMS99]  of  the  partial  differential 
equation  with  a  certain  range  of  initial  data,  in  the  limit  that  the  two  wave  speeds 
coincide.) 

Figure  6  shows  stability  curves  (graphs  of  A  =  A (k))  for  three  of  the  compressive 
waves  from  U-  to  u+  (the  Erst  three  in  the  ordering  of  [BMS99]  that  are  stable  to  one¬ 
dimensional  perturbations),  for  the  undercompressive  wave  and  for  the  trailing  wave. 
The  results  for  the  first  and  third  traveling  wave  are  shown  by  solid  lines,  with  addi¬ 
tional  symbols  (bullets,  pluses  and  diamonds)  for  the  first,  third,  and  fifth  traveling  wave. 
A  dot-dash  curve  represents  the  eigenvalue  for  the  undercompressive  wave,  and  a  dashed 
curve  for  the  trailing  compressive  wave.  We  do  not  use  a  solid  line  for  the  fifth  compres¬ 
sive  wave  since  it  would  coincide  with  the  undercompressive  and  trailing  compressive 
wave  curves.  From  the  figure,  we  immediately  observe  that  the  compressive  waves  all 
have  a  range  of  k  >  0  for  which  they  are  unstable,  whereas  the  undercompressive  profile 
is  linearly  stable  against  transverse  perturbations. 

Figure  6  also  shows  the  following  interesting  phenomenon.  Let  A {k\  i)  denote  the 
stability  curve  for  the  (2 i  —  l)th  traveling  wave  (recall  that  the  even  numbered  waves  are 
unstable  to  one  dimensional  perturbations),  and  let  \uc(k)  and  Xtr{k)  denote  the  domi¬ 
nant  eigenvalues  for  the  undercompressive  and  trailing  compressive  waves,  respectively. 
In  Fig.  6,  observe  that  the  graphs  of  A(/c;  3)  and  of  ma,x{Xuc(k),  Xtr(k)}  are  virtually 
indistinguishable.  This  observation  suggests  that 

lim  X(k ;  i)  —  max{Auc(A;),  Xtr(k)}.  (3.28) 

i— >oo 

This  unusual  limiting  behavior  may  be  understood  as  follows:  for  large  i,  the  (2i  —  l)th 
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Figure  6:  The  stability  curves  for  general  k.  The  results  for  the  three  simple  compressive 
waves  are  given  by  bullets,  plus  and  diamonds.  The  curves  for  waves  one  and  three  are 
filled  in  by  solid  lines.  For  the  trailing  compressive  wave,  the  stability  curve  is  given  by 
a  dashed  line,  and  a  dot-dashed  line  for  the  undercompressive  base  profile. 


traveling  wave  has  a  well-separated  leading  and  trailing  front  which  closely  resembles  a 
double  shock  composed  of  a  leading  undercompressive  shock,  and  a  trailing  compressive 
shock  having  the  same  speed,  each  of  which  have  traveling  wave  profiles.  Correspond¬ 
ingly,  the  compressive  traveling  wave  trajectories  in  phase  space  (i.e. ,  trajectories  from 
the  middle  equilibrium  M  to  the  bottom  equilibrium  B)  are  approximated  by  the  union 
of  a  trajectory  from  M  to  T  of  the  trailing  compressive  wave  and  the  trajectory  from  T 
to  B  of  the  undercompressive  trajectories  at  these  parameter  values.  Hence  the  eigen¬ 
values  for  these  approximately  composite  waves  are  approximately  given  by  the  union 
of  the  eigenvalues  of  the  component  undercompressive  and  trailing  waves.  These  gen¬ 
eral  principles  are  familiar  from  the  study  of  multi-hump  solitary  wave  solutions  in  the 
reaction-diffusion  and  nonlinear  optics  literature,  [AJ1,  AGJ1,  L,  S.l,  S.2]. 

Though  the  proofs  (from  the  solitary  wave  theory)  in  general  break  down  for  systems 
of  conservation  laws,  due  to  the  lack  of  a  spectral  gap  at  A  =  0,  in  the  present,  scalar 
setting  they  can  be  applied  unchanged,  after  the  usual  weighting  transformation  to  re¬ 
cover  a  spectral  gap.  (Alternatively,  one  could  recall  from  Section  2.2  the  fact  that  <S+ 
and  U~  remain  spectrally  separated  at  k  —  0,  A  =  0,  and  carry  out  a  direct  proof  using 
Evans  function  methods.  These  issues  are  discussed  in  [Z.3]). 

Likewise,  the  eigenfunctions  tend  to  be  superpositions  of  the  eigenfunctions  (shifted 
so  they  essentially  do  not  overlap,  in  keeping  with  the  trajectories  themselves)  of  the 
undercompressive  and  trailing  compressive  waves.  The  growth  rate  of  the  composition 
is  then  governed  by  the  maximum  of  the  growth  rates  of  each  part. 

This  observation  suggests  the  following  property: 
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At  k  —  ki,  where  the  two  stability  curves  \uc(k),  Xtr(k)  intersect,  we  might  expect 
a  loss  of  smoothness  in  A (k]i)  for  i  — >  oo,  and  possibly  for  large  enough  finite  i.  The 
behavior  near  k  —  k\  is  explored  further  in  3.4.3. 

3.4.2  Comparison  with  long-wave  asymptotics 

Next,  we  examine  the  range  0  <  k  <  0.1  in  greater  detail,  and  compare  A (k)  for  each 
compressive  wave  with  the  asymptotic  form  A1  A;2,  the  coefficients  A1  being  calculated 
separately  for  each  wave  from  the  formula  (3.20).  The  results  are  shown  in  the  log-log 
plot  of  Figure  7  (left).  In  this  figure,  symbols  correspond  to  computed  values  of  X(k) 
with  triangles  for  the  trailing  wave,  and  circles,  pluses,  and  diamonds  for  waves  one  three 
and  five  respectively.  The  times  symbols  on  the  right  are  for  the  undercompressive  wave. 
The  solid  lines  represent  A1^2  for  compressive  waves  one  and  three.  The  line  for  wave 
five  is  indistiguishable  from  that  for  wave  three.  The  dashed  line  represents  A1^2  for  the 
trailing  wave  and  the  dot-dashed  line  for  the  undercompressive  wave. 

In  Figure  7  (left),  we  note  that  the  graphs  of  X(k]i),i  =  1,2  clearly  approach  the 
asymptotic  curves  as  k  — >  0,  and  stay  reasonably  close  throughout  the  range  0  <  k  <  0.1. 
However,  the  graph  of  X(k,  3)  switches  over  in  this  range  of  k  from  the  asymptotic  curve 
to  the  curve  for  the  compressive  trailing  wave.  We  conclude  the  following  property: 

Fig.  7  (left)  suggests  that  the  coefficient  A)  for  the  (2 i  —  l)th  compressive  traveling 
wave  has  a  limit  and  that 

lim  Xi  ^  Xtr. 

% — ¥  CO 

Moreover  the  numerics  suggest  that  for  larger  i,  the  leading  order  long-wave  asymptotics 
agrees  with  the  spectrum  for  a  smaller  range  of  k  near  zero.  We  conjecture  that  the 
analytic  expansion  breaks  down  in  the  limit  as  k  — >  oo  where  we  are  essentially  linearizing 
about  a  composite  wave. 

For  the  undercompressive  wave,  the  linear  stability  analysis  for  general  k  can  be 
done  with  the  same  numerical  method  as  for  the  compressive  case.  To  compute  the 
coefficient  of  the  long  wave  expansion  A1,  we  have  to  start  with  (3.18),  and  use  the  non¬ 
constant  eigenfunction  i r°  satisfying  (3.15)  with  boundary  conditions  (3.26).  Note  that 
the  argument  presumes  that  A1  <  0,  which  can,  in  part,  be  justified  here  by  referring  to 
the  numerical  results  for  general  k. 

The  adjoint  problem  (3.15),  (3.26)  was  solved  by  discretizing  the  equation  with  finite 
differences  similar  as  for  the  general-A;-problem,  then  solving  for  7T°  iteratively  in  a  process 
similar  to  the  iterative  vector-iteration  process  used  to  calculate  eigenvectors  for  known 
eigenvalues  [Ke,  PTVF].  Specifically,  let  p  —  (pi),  i  —  1...N  denote  the  vector  of 
grid-values  approximating  7r°,  and  L^isc  the  discrete  N  x  N- matrix  representation  of  the 
continuous  operator  given  by  (3.15),  and  the  discrete  boundary  conditions 

P0=P2~  2pi  +  Po  -  0,  pn  -  pn- 1  =  Pa  -  2pn-i  +  pn—2  =  0.  (3.29) 

Then,  successive  approximations  to  7T°  are  obtained  by  solving 

_  discP°ld 
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Figure  7:  (left)  Stability  curves  for  the  compressive  profiles  for  small  k  in  a  log-log  plot. 
The  symbols  (circles,  pluses,  diamonds)  represent  the  numerically  calculated  A (k)  (for 
the  first,  third  and  fifth  wave,  resp.).  Solid  lines  show  the  A1  A;2  for  waves  one,  three  and 
five.  The  dashed  line  shows  A1  A2  for  the  trailing  wave,  (right)  Stability  curve  for  the 
undercompressive  shock.  Symbols  denote  the  values  computed  for  A (k)  using  the  code 
for  general  k,  the  dot-dashed  line  represents  the  long-wave  approximation. 


where  j |  —  | |oo  denotes  the  discrete  maximum  norm.  Note  that,  even  though  Lq  is  singular, 
the  spectrum  of  the  discretized  operator  will  typically  not  include  the  eigenvalue  zero, 
i.  e.  this  operator  is  invertible  and  the  above  step  can  be  carried  out.  After  very  few 
iterations,  in  fact  we  found  in  our  numerical  trials  after  the  first  iteration,  the  sequence  of 
discrete  p  becomes  stationary.  The  vector  thus  obtained  is  eigenvector  for  the  eigenvalue 
with  the  smallest  modulus,  and  hence  a  good  approximation  of  7T°.  For  the  special 
situation  of  equation  (3.27),  the  resulting  7T°(x)  is  shown  in  Figure  8.  Also  note  that 
the  normalization  of  p  —  0  with  respect  to  the  maximum  norm  at  each  step  prevents  it 
from  converging  to  the  zero  vector;  on  the  other  hand,  the  discrete  boundary  conditions 
(3.29)  exclude  convergence  to  a  non-zero  constant  solution. 

This  7T°  is  then  used  to  evaluate  (3.18).  The  resulting  plot  for  A1^2  is  shown  in 
Figure  7  (right)  as  a  dot-dashed  line.  The  agreement  with  the  numerical  values  A (k)  is 
excellent . 

3.4.3  The  spectrum  for  higher  order  compressive  waves 

We  now  present  some  more  detailed  computations  showing  the  structure  of  the  higher 
order  compressive  waves.  Using  an  iterative  method  similar  to  the  one  used  to  compute 
7T°,  we  can  compute  the  top  two  eigenvalues  of  the  linearized  operator  corresponding  to 
traveling  waves  three  (i  —  2)  and  five  (i  —  3). 

Figure  9  shows  the  top  two  eigenvalues  and  corresponding  eigenfunctions  for  wave 
number  three  (i  —  2).  On  the  left,  the  solid  line  shows  the  top  eigenvalue  and  the 
dotted  line  shows  the  bottom  eigenvalue.  The  long-dashed  line  shows  the  dominant 
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Figure  8:  The  non-constant  left  eigenfunction  for  the  undercompressive  profile. 


Figure  9:  (left)  Spectrum  of  the  first  two  eigenvalues  (shown  as  solid  and  dotted  lines) 
for  the  third  compressive  wave  (i  —  2).  The  dot-dashed  line  and  long-dashed  line  show 
the  dominant  eigenvalue  for  the  undercompressive  wave  and  trailing  wave,  respectively. 
The  inset  is  an  enlarged  view  of  a  region  near  k  =  0,  marked  by  a  box  with  dotted  lines 
in  the  outer  graph,  (right)  The  eigenfunctions  for  the  first  two  eigenvalues  of  traveling 
wave  three  {i  =  2)  at  k  =  0. 
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Figure  10:  (left)  Spectrum  of  first  two  eigenvalues  (shown  as  solid  and  dotted  lines)  for 
the  fifth  (i  —  3)  compressive  wave.  The  dot-dashed  line  and  long-dashed  line  in  the 
inset  show  the  dominant  eigenvalue  for  the  undercompressive  wave  and  trailing  wave, 
respectively,  (right)  The  eigenfunctions  for  the  first  two  eigenvalues  of  traveling  wave 
five  (i  —  3)  at  k  —  0. 


eigenvalue  for  the  trailing  wave  and  the  dot-dashed  line  the  dominant  eigenvalue  for  the 
undercompressive  wave.  It  is  clear  for  traveling  wave  three  that  the  two  eigenvalues  are 
well  separated  and  that  the  second  eigenvalue  approaches  a  small  but  negative  value 
as  k  — >  0.  On  the  right  are  the  eigenfunctions  for  the  top  two  modes  at  k  =  0.  The 
eigenfunction  with  A  =  0  is  shown  as  the  solid  line;  it  agrees  precisely  with  the  translation 
mode  ux  shown  via  circles.  The  eigenfunction  corresponding  to  the  negative  eigenvalue 
is  shown  via  a  dotted  line;  it  corresponds  to  infinitesimal  but  unequal  shifts  of  the 
leading  and  the  trailing  edges  of  the  wave,  indicated  by  the  observation  that  this  wave 
is  approximately  cux  for  c  >  1  near  the  trailing  part  of  the  wave  and  ux  near  the  leading 
part  of  the  wave.  The  small  negative  value  of  this  eigenvalue  at  k  —  0  reflects  the  fact 
that  the  wave  is  just  stable  against  one-dimensional  perturbations  towards  one  of  its 
neighboring  waves. 

Figure  10  shows  the  top  two  eigenvalues  and  corresponding  eigenfunctions  for  wave 
number  five  (i  —  3).  On  the  left,  the  solid  line  shows  the  top  eigenvalue  and  the 
dotted  line  shows  the  bottom  eigenvalue.  The  long-dashed  line  showing  the  dominant 
eigenvalue  for  the  trailing  wave  and  the  dot-dashed  line  the  dominant  eigenvalue  for 
the  undercompressive  wave  are  only  shown  in  the  inset  as  their  agreement  with  the  top 
two  eigenvalues  of  wave  number  five  is  very  good.  Unlike  traveling  wave  three,  the  two 
eigenvalues  are  not  well  separated  and  appear  to  cross  in  the  region  marked  with  a  box 
and  blown  up  in  Figure  11.  However,  as  for  traveling  wave  three,  the  second  eigenvalue 
approaches  a  small  but  negative  value  as  k  — >  0.  The  inset  shows  an  enlarged  view  of  a 
region  near  k  —  0,  delineated,  in  the  outer  graph,  by  a  small  box  with  dotted  lines.  For 
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Figure  11:  (left)  Enlarged  view  of  the  real  part  of  the  Erst  two  eigenvalues  (show  as  solid 
and  dotted  lines)  for  the  fifth  compressive  wave.  The  dot-dashed  line  and  long-dashed 
line  show  the  dominant  eigenvalue  for  the  undercompressive  wave  and  trailing  wave, 
respectively,  (right)  The  imaginary  part  of  the  eigenvalues  show  on  the  left. 


a  blow-up  of  the  region  near  k  —  k\  marked  by  the  larger  box  with  solid  lines,  see  fig.  11. 

On  the  right  are  the  eigenfunctions  for  the  top  two  modes  at  k  —  0.  The  eigenfunction 
with  A  =  0  is  shown  as  the  solid  line;  again  it  agrees  precisely  with  the  translation  mode 
ux  shown  via  circles.  The  eigenfunction  corresponding  to  the  negative  eigenvalue  is 
shown  via  a  dotted  line;  it  again  corresponds  to  infinitesimal  but  unequal  shifts  of  the 
leading  and  the  trailing  edges  of  the  wave. 

Figure  11  shows  an  enlarged  view  of  Re  A (k)  and  Im  A (k)  near  the  point  where  the 
two  eigenvalues  appear  to  cross  in  Figure  10.  The  line  styles  carry  over  from  that  figure. 
As  k  approaches  the  intersection  point,  the  top  and  second  eigenvalues  coalesce  giving 
rise  to  a  complex  conjugate  pair  of  eigenvalues,  which  reseparate  into  real  eigenvalues 
after  the  intersection  point.  The  imaginary  part  of  the  eigenvalues  is  shown  to  grow  from 
zero  and  then  shrink  back  to  zero  during  this  transition  (see  right  figure).  The  vertical 
dotted  lines  on  each  graph  have  exactly  the  same  k  values  and  delineate  the  bifurcation 
points  in  k. 
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